INCORPORATING STRUCTURAL PRIOR INFORMATION AND SPARSITY INTO EIT USING PARALLEL LEVEL SETS

. EIT is a non-linear ill-posed inverse problem which requires so-phisticated regularisation techniques to achieve good results. In this paper we consider the use of structural information in the form of edge directions coming from an auxiliary image of the same object being reconstructed. In order to allow for cases where the auxiliary image does not provide complete information we consider in addition a sparsity regularization for the edges appearing in the EIT image. The combination of these approaches is conveniently described through the parallel level sets approach. We present an overview of previous methods for structural regularisation and then provide a variational setting for our approach and explain the numerical implementation. We present results on simulations and experimental data for diﬀerent cases with accurate and inaccurate prior information. The results demonstrate that the structural prior information improves the reconstruction accuracy, even in cases when there is reasonable uncertainty in the prior about the location of the edges or only partial edge information is available.


(Communicated by Samuli Siltanen)
Abstract. EIT is a non-linear ill-posed inverse problem which requires sophisticated regularisation techniques to achieve good results. In this paper we consider the use of structural information in the form of edge directions coming from an auxiliary image of the same object being reconstructed. In order to allow for cases where the auxiliary image does not provide complete information we consider in addition a sparsity regularization for the edges appearing in the EIT image. The combination of these approaches is conveniently described through the parallel level sets approach. We present an overview of previous methods for structural regularisation and then provide a variational setting for our approach and explain the numerical implementation. We present results on simulations and experimental data for different cases with accurate and inaccurate prior information. The results demonstrate that the structural prior information improves the reconstruction accuracy, even in cases when there is reasonable uncertainty in the prior about the location of the edges or only partial edge information is available.
1. Introduction. Consider the forward problem (1) y meas = A(f ) + n , where y meas ∈ Y = R M is a discrete set of observations, f ∈ X is an image in a solution space X, A : X → Y is the forward operator mapping from solution space to data space and n ∈ Y is measurement noise. Image reconstruction consists of recovering an estimate of f from y meas , i.e. in solving the inverse problem (1). Within this general class of problems we may identify image processing problems such as denoising, deblurring, in-painting and super-resolution, wherein the spaces X ⊇ Y are discrete sets of pixels, and tomographic problems wherein X is a function space and A is a continuous-discrete operator that may be weakly or strongly illposed, and either linear or non-linear.
A general approach to such inverse problems is a regularised output-least-squares scheme whereby the estimate is the solution to an optimisation problem (2) f * = arg min with Γ being the covariance matrix of the measurement noise and Ψ(f ) a regularisation functional that may be formally considered as the negative logarithm of a prior probability density model for f , and α > 0 the regularisation parameter. Note that the form of (2) assumes that the noise n is modelled as zero mean Gaussian n ∼ N (0, Γ). Many different forms of regularisation have been considered in the literature. A topic of high interest is the use of priors based on sparsity assumptions; in particular priors that promote images with a sparse description of their edges have played a significant role in many areas of image reconstruction. Separately, there exists a class of priors assuming the availability of an auxiliary image p and where some assumptions are made about the correlation of the two images, for example in terms of similar structural and/or statistical information. Recently, we reviewed some strategies for the combination of these approaches to define a sparsity promoting prior that takes into account the edge information of both the image f that is being reconstructed, and the auxiliary reference image p [18]. In that paper we applied this idea to the problem of image reconstruction for Positron Emission Tomography (PET) using a co-registered anatomical image, acquired with Magnetic Resonance Imaging (MRI), as the auxiliary image and a novel idea for an edge-promoting prior using a parallel level sets strategy where structural similarity of two images is measured by differences in the alignment of the level sets of the images. The core of the method is an anisotropic regularization functional Ψ(f ) which promotes alignment of the level sets of the unknown image f with the level sets of the reference image p. In particular, Ψ(f ) is designed so that it poses small penalty for changes in f perpendicular to level sets where the reference image p exhibits large changes (e.g., organ boundaries). In this paper we investigate the application of this prior to Electrical Impedance Tomography (EIT), a non-linear, severely ill-posed reconstruction problem. Whereas in this paper our motivation is to investigate both the feasibility and usefullness of incorporating structural information as a particular regularisation scheme, we note that, in general, the concept of multi-modality imaging is of increasing interest, because different modalities provide complementary information regarding the patient being studied. In particular, the availability of e.g. an MRI or CT image does not obviate the usefullness of also having an EIT image, because the latter provides information on physical parameters (specifically, electrical conductivity) which the standard modalities cannot provide.
In EIT the problem is to reconstruct an interior conductivity from knowledge of full or partial Cauchy data on the boundary (i.e. the Dirichlet-to-Neumann map). In general this problem is known to possess a logarithmic-type stability [1,35] since high spatial-frequency components in the interior are increasingly strongly damped in the forward mapping from interior conductivity to boundary voltage measurements. It is this fact that gives rise to the intrinsically low-resolution of reconstructed EIT images. The role of regularisation is to introduce prior knowledge on the distribution of possible solutions which is formally stated in the Bayesian framework. In particular, a prior which emphasises the conductivity image to have well defined edges whose location is consistent with a known or partially known geometric structure is the main aim of this work. This is motivated by ongoing efforts for development of EIT based bedside monitoring of brain stroke and lung function in emergency care [21,2,16,34]. In both of these applications, the patient specific anatomy (head or chest) is usually available from a CT scan that has been taken for diagnosis before admission to the emergency care where the bedside monitoring takes place, and our hypothesis is that utilizing the CT information in the EIT reconstruction in a suitable way could yield significant improvement in the accuracy of the EIT based imaging.
In practice, edge information might be obtained for example using edge-detection algorithms on the auxiliary image p. Since in general such algorithms may be subject to error, it is important to consider the effect of inaccurate edge location information on the performance of our methods. Also, in cases when there is uncertainty about the location of the boundaries, one could potentially utilize a spatial weighting model in the parallel level set regularization that accounts for uncertainty in the locations of the edges in the prior. Such a weighting model could be constructed, for example, using anatomical atlases when patient specific CT is not available or deliberate blurring of the available anatomical reference image to account for motion of the boundaries. We will consider simple examples of incorrect edge location and uncertain edge information in the numerical examples.
1.1. Related work. In this paper, we investigate utilization of structural prior into EIT using the parallel level set strategy where structural similarity is measured by differences in the alignment of the level sets of the images. However, structural priors can be formulated in variety of ways and find use in many settings. We will briefly review other application domains and options for how to include a-priori structural information.
1.1.1. Application domains. Over the years there have been many different areas where structural regularisation has been applied. The dominant domain for structural priors has been medical imaging. Since the early 1990's [32] high resolution anatomical knowledge-obtained either from MRI or computer assisted tomography (CT)-has been incorporated into the reconstruction of a low resolution functional image, in particular in PET and Single Photon Emission Computed Tomography (SPECT), and it still continues to be an active area of research, e.g. [10,43,25,4,18]. This application domain has gained more interest recently with the advent of combined PET-MRI [15]. Next to nuclear medicine, similar techniques have been used to exploit the predicted correlation between different tissue contrasts in MRI, e.g. [17] and in functional MRI and contrast enhanced dynamic MRI [38]. In [6] MRI is used for construction of structurally guided regularization for magnetic particle imaging (MPI). In [2] structural prior information is exploited in the EIT Dbar reconstruction method. Beside the dominant field of medical imaging we would also like to point to two other areas which have been often overlooked in terms of structural priors: remote sensing and image processing.
Remote sensing often deals with images that are of high spectral resolution but have limited spatial resolution. When a panchromatic image with high spatial resolution is available, these two images can be fused to create an image of high spectral and high spatial resolution. In remote sensing this technique is often referred to as pansharpening and can be achieved with structural priors, e.g. [5,36,20,12] and references therein.
Structural priors have also been used in image processing despite the fact that there is usually no a priori structural information available. Instead, in an iterative process, the structure of an image is estimated and the image is restored. This technique has been applied to estimate a spatially varying regularization parameter field, e.g. for total variation [40,22,29,33] and total generalized variation [11], or anisotropic structures [23,33,19].
All these application fields have in common that the forward operator is modeled linearly; in contrast, in the EIT setting of this paper we are facing a non-linear inverse problem.
1.1.2. Types of structure. Structural priors come in a variety of forms. Most of the earlier works on structural priors are isotropic, meaning that only the locationrather than the direction-of structure has been taken into account. This is mostly achieved by spatially varying weights, e.g. [40,3,22,29,33,17,11], but can also be achieved in a joint/vectorial total variation fashion [18].
More information can be extracted by anisotropic structures. While the isotropic priors typically depend only on the local magnitude of the gradient of the prior image, the anisotropic models also take into account the alignment of the edges. The anisotropy can be defined globally [7] or spatially varying [27,5,28,23,33,17,18,19,12] as in this contribution. A common way of modelling anisotropy is by a spatially varying vector-or matrix-field that defines preferred directions. Positive correlation of structure can be incorporated with Bregman distances [5,28] and general correlation with projections [27,5,23,33,17,18,19,12]. Anisotropic priors are found to be more effective in applications, see e.g. [17], thus we solely focus on these in this study.
Some of the structural priors are defined in a discrete setting as Markov random fields where it is difficult to make the distinction between isotropic and anisotropic structures [32,10].
1.1.3. Algorithms for structural regularisation. Most of the aforementioned priors are generalizations of total variation and as such can make use of all the algorithms that handle normal total variation. For instance, gradient based algorithms can be used if the problem has been smoothed, e.g. [23,18], fast gradient projection on the dual problem [8] has been used in [17] and primal-dual algorithms [13] in [33,19]. As state-of-the-art primal-dual or proximal methods are not directly applicable to our non-linear (and thus non-convex) setting, we focus here on gradient based algorithms for smooth objectives and use smooth approximations of non-smooth functionals.

1.2.
Organization of the paper. The rest of the paper is organised as follows: in section 2 we derive a notation for variational regularisation based on the gradient of the unknown with structural information, and we derive the Euler-Lagrange operator in terms of a second order PDE of diffusion type. This allows us to relate regularisation to various well-known diffusive flows as used in image enhancement methods. In section 3 we discuss the EIT problem and explain the numerical implementation of the structured regularisation for EIT using the parallel level sets approach. In section 4 we present results from simulations and experimental data from a laboratory setup. In section 5 we give conclusions.
2. Variational regularisation. A classical way to approach the optimisation problem (2) is through variational methods applied to the functional Φ(f ) whose first variation defines a direction of steepest descent Since Φ is defined as the sum of a data fit term and a penalty term, we can consider their first variation forms separately, which leads to the concept of variational regularisation. Therefore in this section we first review some existing variational regularisation schemes and their corresponding derivatives in section 2.1, before considering how they are modified in the presence of structural prior information in section 2.2. Finally in section 2.3 we combine these with the data fit term to obtain an overall nonlinear reconstruction scheme.
2.1. Image regularisation and diffusive flow. We restrict our discussion to regularisation functionals of the form wherein a local mapping ψ : R → R modifies the weighting of the local image gradient ∇f . The corresponding first variation leads to an Euler-Lagrange equation that defines a PDE of diffusion type for image flow where the local diffusivity function κ = ψ ( ∇f )/ ∇f controls the flow of intensity between pixels. This can equivalently be stated in terms of the Euler-Lagrange operator L as defined in (5), Indeed different choices of the mapping ψ give rise to different scale spaces f (x, t) under the evolution (7). Some commonly used choices for ψ, with the resulting form of κ, are given in table 1. Note that, apart from the (non-smooth) total variation function, these priors are chosen so that κ ∈ [0, 1] which simplifies the analysis of the stability of (7). Furthermore in table 1 we introduced the threshold value T to explain the smoothing of the mapping ψ in terms of the dimensionless quantity ∇f /T : T indicates the transition from edges arising from noise, and those from true structure in the image, with the expectation that κ → 1 below the threshold and (7) tends to an isotropic blurring with a local Gaussian kernel.

Incorporating structural information with parallel level sets.
A seminal paper demonstrating the use of structural prior information by "parallel level sets" was [27]. In that paper the 1 st -order Tikhonov scheme for Ψ was used (i.e. using ψ(t) = 1 2 t 2 ) in (4) with a fixed (symmetric) tensor field B(p) (or B for short), incorporating structure from an auxiliary reference image p, used to modify the metric of the local norm In a nutshell, the tensor B can be chosen such that structures (represented by level sets) in f and p are favoured to be aligned (or parallel), thus coining this method "parallel level sets". Such a tensor B is for instance given in terms of the local Table 1. Examples of ψ for different regularisation schemes in variational form. ψ is the mapping which defines the penalty for the gradient magnitude in (4) and κ is the corresponding local diffusivity function in (6).
gauge coordinatesν ⊥ normal, andν tangent to the level sets of p, with B diagonal in the rotated system Here Λ is an anisotropy matrix controlling the relative flow normal and tangent to the level sets of p and 0 ≤ γ ≤ 1 is an edge indicator function tending to zero as ∇p → ∞. The diffusive flow resulting from the first variation of (8) is T is the gradient operator in the local gauge coordinates. Note that this flow reduces to isotropic heat flow if p is locally flat, i.e. γ = 1.
With this notation we can write with L = RΛ 1/2 and (8) can also be written To combine this approach with the sparsity concept, we generalise (11) to whose Euler-Lagrange operator can be defined as As we will see in the results section, if ψ is chosen to be TV-like, this scheme reduces to total variation (or a smoothed approximation thereof) for the regularisation of f when p is flat. In the results section, we consider structured regularization using the smoothed TV functional (see table 1

Regularised reconstruction.
We are now ready to introduce the structured regularisation into the problem (2). Putting the form in (8) into (2) and using a numerical implementation of the lagged Gauss-Newton algorithm [42], leads to an iterative scheme of the form: where A n : X → Y is the Fréchet derivative of the forward operator at f (n) , (A n ) * : Y → X its adjoint, L n = L(f n ) the Euler-Lagrange operator at f (n) and s (n) is a step length parameter obtained from a line search algorithm which finds an optimal step length along the current search direction by solving a 1D minization problem [9,37].
3. EIT. In the following we take the particular notation f → σ to indicate that we consider the image to be a map of electrical conductivity.

Forward model of EIT.
Let Ω be a bounded domain in R d (d = 2, 3) with smooth boundary ∂Ω. In an EIT experiment, a set of N el contact electrodes are attached on the boundary ∂Ω of Ω. A set of electric currents, called current patterns, are injected via the electrodes into the body Ω and corresponding voltages are measured using the same electrodes. We model these measurements with the complete electrode model [14,39] ∇ · (σ(r)∇u(r)) = 0 , r ∈ Ω (17) where r ∈ Ω is the spatial coordinate, σ(r) is the conductivity, u(r) is the electric potential inside Ω, U and I are the potential and current at electrode e , respectively, z is the contact impedance between the electrode e and the body Ω, and ν ∂Ω denotes the outward unit normal vector on the boundary ∂Ω. The charge conservation law and the fixing of the zero potential are implied by The existence and uniqueness of the solution of the model (17)(18)(19)(20)(21) was proven and its variational form derived in [39].
We use a finite element method (FEM) for solving numerically the forward problem (17)(18)(19)(20)(21). In the FEM approximation, the domain Ω is divided into M e disjoint elements joined at M n vertex nodes. The potential u and electrode potentials U ∈ R N el satisfying the variational form [39] of (17)(18)(19)(20)(21) are approximated as where the functions φ i are the piece-wise linear nodal basis functions of the finite element mesh and vectors w j ∈ R N el are chosen such that condition (21) holds.
In this study, the conductivity σ(r) is approximated as a piecewise constant function on a grid of regular square pixels Ω ⊂ N k=1 G k , leading to an approximation with N pixels within or partially intersecting Ω and where χ is the characteristic function taking value 1 on all points in the set of its argument. The pixel values of the conductivity are concatenated into a vector σ = (σ 1 , . . . , σ N ) T ∈ R N . Using these approximations and the FEM implementation described in [41], the numerical forward solution for each current injection is obtained by solving a (M n + N el − 1) × (M n + N el − 1) system of equations. Let V ( ) denote the voltage measurements corresponding to the th current pattern I ( ) and U ( ) (σ) denote the respective FEM based forward solution. Let the EIT experiment consist of K current injections. We denote the measurement vector and forward model for the whole EIT experiment by y meas = (V (1) , . . . , V (K) ) T ∈ R KN el and A(σ) = (U (1) (σ), . . . , U (K) (σ)) T ∈ R KN el . The Fréchet derivative of A(σ) is found by solving the adjoint problem to (17)(18)(19)(20)(21) for u * and forming the inner product ∇u * j , G k ∇u i for j = 1, . . . , N el , i = 1, . . . , K, k = 1, . . . , N , for more details see [26].
The measurement noise in EIT experiments is commonly modeled as additive noise as in (1) with block-diagonal covariance matrix Γ.
The minimization was carried out using an implementation of the Lagged Gauss-Newton optimization method (16), where the line search was implemented using a bounded minimization algorithm such that positivity of the conductivity is enforced. Two different choices of the regularization functional were considered; SH 1 structured smoothness regularization [27], Ψ(σ) given by (11). STV structured (smoothed) TV regularization, Ψ(σ) given by (15).
In both cases, the tensor field B(p) was of the form (10), whereν ⊥ and γ were selected asν with the parameter T p > 0. We want to remark that these choices lead to a continuous tensor field B(p), in the sense that lim ∇p→0 B(p) = I.
In the discretization of the regularization functional, the gradient of the conductivity was approximated using a forward differencing scheme ∇σ k = l −1 (σ h − σ k , σ v − σ k ) T where l is the distance between the center points of the neighboring pixels and σ h and σ v are the values of the conductivity in the neighbor pixels of pixel k in the (forward) horizontal and vertical directions. The gradient of the reference image p was approximated by forward differences as well.
In all cases the original reference image p was in finer resolution than the N -pixel representation (24) used for the conductivity. Thus, the reference image was mapped to the same piecewise constant pixel basis that is used for σ before evaluation of ∇p.

4.
Results. In this section, the structured regularisation is evaluated using simulated data and experimental EIT data from a laboratory setup.
In each test case, reconstructions with conventional (no structure) regularisation were computed as a reference for the structured regularisation. This was obtained by setting p ≡ const., leading to B = I. The reconstructions with the structured regularisation were computed using both correct edge information and partial edge information. The correct edge information refers to a case where the reference image contains correct edge information in the sense that ∇σ true and ∇p are parallel everywhere in the domain Ω, but the signs and/or values of ∇σ true and ∇p may be different. The partial edge information refers to a case where the edges in σ true and p are partially different, or a case with uncertainty about the location of the boundaries. Specially, the partial edge information that is used in the simulated test cases has the following features: (Case 1). Both the conductivity σ and the reference image p have one inclusion in a location where the other image is constant. This could be considered as an example of a case with one organ that has contrast only in conductivity and one organ that has contrast only in the reference image. This example is important to demonstrate that the structural prior does not promote false information. (Case 2). The inclusions in the partial edge information are in the correct location and roughly of correct size but their topology is incorrect. This case can be considered as an example of where the resolution of the auxiliary image p is at a different scale to the conductivity σ.
(Case 3). The partial edge information model simulates a case with uncertainty on the locations of the edges. The uncertainty in the locations is modelled by using a piecewise constant weighting function γ which gives equal weighting over the area where the boundaries are expected to be. The results include a study on how the magnitude of the uncertainty in the locations affect the reconstruction quality. (Case 4). Similar to Case 3 with the exception that the true conductivity contains one nested inclusion which is not present in the reference image p. The case can be viewed as a simplistic example of imaging the lungs in a case where one of the lungs contain an anomaly that is not present in the prior. The regularisation parameter α was selected experimentally. For simulated test cases the selection was carried out by computing the mean of the relative L 2 error where P is an interpolation mapping which maps the reconstructed conductivity to the same grid as the true conductivity, with different values of α over 20 realizations of data y meas for the first test case in figure 1, and then selecting α which gave the smallest reconstruction error for the conventional (no structure) case p = const. This procedure resulted to values of α = 5 for the SH 1 regularisation and α = 0.5 for the STV regularisation. The same fixed values of α were then used in all the simulations with all the different choices of reference image p. This choice was made so that the results for the conventional (no structure) regularisation, which is shown as reference for the structured regularisation, would be optimal. The values of α for the experimental test case were selected manually based on the conventional (no structure) regularisation, and then the same fixed values were used with all the different choices of p. For the SH 1 regularisation the value was α = 25 and for the STV regularisation α = 2. 4.1. Simulated data. The domain Ω in the simulations was a circle with radius 14 cm with N el = 16 equispaced, 26 mm wide electrodes located on the boundary ∂Ω. For the simulation of the measurement data, the domain Ω was divided into a mesh of 44 006 triangular FEM elements with 22 820 nodes. The number of regular pixels for the representation of the true target conductivity σ true was N = 13 677. The data was simulated using a pairwise current injection paradigm using 15 current patterns such that electrode e 1 was the source in each current pattern and electrode e 1+t (t = 1, . . . , 15) the sink. For each current pattern, voltages were recorded between pairs of neighboring electrodes, leading to 240 voltage measurements for each test case. The electrode contact impedances were set to z ≡ 0.001 in the simulations and were treated as known parameters in the computation of the reconstructions. Zero mean Gaussian white noise with standard deviation σ e = 10 −4 was added to the simulated voltages. With this choice of σ e , the signal-to-noise ratio (SNR) in the simulated voltage measurements were in the range from 68 dB to 70 dB, and the ratio of σ e divided by the absolute values of the noiseless voltages were in the range from 0.01 % to 1.1 %, representing noise statistics that are similar to measurements with practical EIT systems, e.g. the KIT 4 EIT system [31].
In the computation of the reconstructions, the domain Ω was divided into M e = 29 932 triangular FEM elements with M n = 15 207 nodes for approximation of u and N = 6 357 rectangular pixels for the representation of the unknown conductivity. Thus, the unknown in the inverse problem was σ ∈ R 6 357 . The results for the first test case are shown in figure 1. The image on the top in the first column shows the true target conductivity σ true . The second column shows results using conventional (no structure) regularisation (p ≡ 1), the third column shows results using structured regularisation using correct edge information and the fourth column results using partial edge information, where the reference image p contains two inclusions, one correctly located and one incorrectly located. The first row in the figure shows the reference images p and the second row the respective weighting functions γ, given by eq. (27) with T 2 p = 0.1. The third row shows the reconstructions using the SH 1 regularisation and the fourth row the reconstructions using the STV regularisation. The reconstruction errors are tabulated in table 2. Figure 2 shows the standard deviations of the reconstructed conductivity with respect the bias of the reconstructed conductivity with different values of the regularisation parameter α for the test case in figure 1. The left image shows the curves for the mean conductivity in the area of the true inclusion in the top of the target conductivity σ true and the right image for the area of the inclusion on the bottom right in the target conductivity. The biases and standard deviations were estimated by computing the reconstructions over an ensemble of 20 independent realizations of the measurement data. The value of α varied in the range [0. 1 20] for the SH 1 regularisation and in the range [0.01 10] for the STV regularisation. The triangle denotes the point corresponding to the smallest value of α in the curves.
As can be seen from figures 1, 2 and table 2, the structured regularisation improves the accuracy over the conventional (no structure) regularisation for both smoothness and TV regularisation. When the correct edge information is used, the reconstructions with both regularisation functionals are very similar, but the STV regularisation yields slightly smaller reconstruction error than SH 1 regularisation. When no edge information or partial edge information is used, the STV regularisation leads to more accurate reconstruction than the SH 1 regularisation-this however can be expected in a case of a piece wise regular target conductivity such as the case considered here. While the inclusion which is on the bottom left in the reference image p but not in σ true , can be seen as a very weak shadow in the SH 1 reconstruction using the partial edge information, it is not present in the STV reconstruction using the partial edge information. Also, another notable difference between the SH 1 and STV regularisation in the partial edge information case is that the STV leads to significantly better recovery of inclusion on the bottom right which is missing from the reference image p. The biases and standard deviations in figure 2 indicate that the structured regularisation leads to more accurate reconstruction of the case in figure 1 than conventional regularisation with any σ true p(r) (no structure) p(r) (correct) p(r) (partial) γ(r) SH1 STV Figure 1. Numerical experiment (Case 1). Rows top to bottom: σ true and reference images p(r) (top row), weighting function γ(r) (second row), reconstructions using SH 1 regularisation (third row) and STV regularisation (fourth row).
value of α. When starting from the biases with the smallest α, the biases decrease in the beginning until reaching a range of suitable values of α. This behavior can be attributed to the fact that the initial α values are too small (under regularisation), leading to a situation where the modelling (discretization) errors cause artifacts to the reconstructed conductivity. When going to (too) large values of α, the biases start increasing (too much regularisation). This behavior is especially evident for the conventional (no structure) regularisation (blue line) where the estimates approach asymptotically zero as the value of α increases. For the structured regularisation, the bias increases more slowly due to different asymptotical behavior which can be explained, loosely speaking, by the different structure of the subset where Ψ(f ) is (essentially) zero. However, the structured regularisation has, for all values of α, smaller than or equal error compared to the conventional regularisation, both for correct and partial edge information. Figure 3 shows results for the second test case. All the estimates were computed similarly as in the first test case. The reconstruction errors are tabulated SD BIAS Figure 2. Plot of standard deviation with respect to bias of the reconstructed conductivity with different values of the regularisation parameter α for the simulation in figure 1. The left image shows the curves for the mean conductivity in the area of the true inclusion on the top in σ true , and the image on the right for the inclusion on the bottom right. The triangle denotes the point corresponding to the smallest value of α in the curves.
in Table 2. The first column shows the true target conductivity, which contains one large inclusion on the left and two side by side circular inclusions on the right. The reconstructions with the conventional (no structure) regularisation are shown in the second column. A notable feature in the conventional reconstructions is that the side by side circular inclusions on the right are reconstructed as one large inclusion with both the SH 1 and STV models. This indicates that the side by side circular inclusions cannot be resolved with the given EIT measurement when using conventional SH 1 or STV regularisation models. The reconstructions with structured regularisation using correct edge information are shown in the third column. Similarly as in the first test case, the structured regularization with correct edge information improves the reconstructions significantly compared to the conventional models. The structured SH 1 and STV reconstructions are visually similar but the STV yields smaller reconstruction error than the SH 1 regularisation. The fourth column shows the results using incorrect edge information. While in the first test case the partial edge information contained one correctly and one incorrectly located inclusion, here the situation is somewhat more challening; the prior shapes are in correct location and of roughly correct size but their topology is incorrect. As can be seen, the incorrect edge information leads here to distortion of the reconstructed conductivity towards the incorrect prior shapes. However, the reconstructions show an increase of the conductivity near the interface of the two circular prior shapes on the left in the area where the true inclusion and the prior shape are different. Similarly, a decrease of conductivity is seen on the right in the area where the true circular inclusions and the prior shape are different. These changes reflect local inconsistency of the prior shapes and they are more evident for the STV regularisation than the SH 1 model. Table 2 shows that the reconstruction errors are smaller σ true p(r) (no structure) p(r) (correct) p(r) (partial) γ(r) SH1 STV Figure 3. Numerical experiment (Case 2). Rows top to bottom: σ true and reference images p(r) (top row), weighting function γ(r) (second row), reconstructions using SH 1 regularisation (third row) and STV regularisation (fourth row).
in both cases with incorrect edge information compared to the conventional (no structure) regularisation. These results suggest that the structured regularisation with incorrect edge information can improve over the conventional model also in cases where the prior shape errors are at the level of measurement resolution power. Figure 4 shows reconstructions for the third test case. The reconstructions with the conventional (no structure) regularisation (p ≡ 1) and the structured regularisation with correct edge information were computed similarly as in the first and second test case. However, the partial edge information was here constructed to simulate a case with uncertainty about the location of the edges. In practice, such a case could occur, for example, when the reference image p would be based on a probabilistic boundary model derived from an anatomical atlas, or the edge location information would be deliberately blurred to account for uncertainty due to movement of the organs, such as the heart and lungs. Here the partial information was constructed by using a smooth reference image (top right in figure 4), which has an approximately linear contrast change over an uncertainty strip where the boundary is expected to be. The smooth reference image was obtained by blurring the correct reference image using an averaging disk filter. To make the weighting of an edge equivalent all over the uncertainty strip, the weighting function in the partial edge information was, instead of using (27), defined as (29) γ(r) = 10 −3 when ∇p(r) > 1 otherwise where is an edge threshold parameter, which in this study was selected manually with value = 0.03. Figure 5 shows the reconstructions with the partial edge information with increasing degree of uncertainty about the location of the edges. The smoothing of the reference image increases to the right, which is seen as widening of the uncertainty strip (in this case, simply the part of domain where ∇p > ). The reconstruction errors for the SH 1 regularisation vary from 8.1 % to 9.8 % from the left to the right in figure 5 (10.6 % for the conventional regularisation using flat p). For the STV regularisation, the reconstruction error changes from 6.3 % to 7.3 % from the left to the right in figure 5 (8.0 % for the conventional regularisation). Figure 6 shows results for the fourth test case. The reconstruction errors are tabulated in Table 2. The test case is otherwise the same as Case 3 in Figure 4 with the exception that the true target conductivity σ true contains a high conducting inclusion inside the larger, low conducting inclusion on the right. The reference images p(r) and weighting functions γ(r) are exactly the same as the ones used in figure 4. This simulation can be considered as a simplistic example of imaging of the lungs where the right lung has an conducting anomaly that is not present in the prior image p(r). The reconstructed images show that the anomaly inside the low conducting inclusion is detected more clearly when structured regularization is employed, especially when using the STV regularization.
The results in figures 4-6 suggest that the structured regularisation can improve the accuracy over conventional (no structure) regularisation even when there is quite large uncertainty in the prior about the exact location or shape of the edges.
The reconstruction times with the structured prior were slightly shorter compared to conventional (no structure) reconstructions for both, the SH 1 and STV, regularizations. As an example, while the reconstruction time for the SH 1 without the structural prior in Case 3 was 1 392 s, the reconstruction time of SH 1 with the correct prior was 1 358 s. The corresponding reconstruction times for the STV regularization were 1 617 s without the prior and 1 330 s with the prior. All these computation times are based on the same implementation of Gauss-Newton algorithm using Matlab on a regular desktop computer, and as such, the computation times are indicative only. However, they do suggest that the use of structured regularization leads to a slight decrease in the computation time compared to the corresponding conventional regularization.

4.2.
Experimental data. The experimental data was measured from a vertically symmetric measurement tank Ω, see figure 7. The cross-sectional shape of the measurement tank was extracted from a CT image of the chest of an adult male. Sixteen equally spaced stainless steel electrodes were attached on the boundary ∂Ω of the tank. The height of the tank Ω was 5 cm. To construct the true conductivity σ, heart and lung shaped inclusions were made of agar and placed in the measurement tank filled with saline. The inclusions were constructed using vertically symmetric moulds. The conductivity of the lung target was approximately 25 % of the conductivity of the saline and the conductivity of the heart target was approximately σ true p(r) (no structure) p(r) (correct) p(r) (partial) γ(r) SH1 STV Figure 4. Numerical experiment (Case 3). Rows top to bottom: σ true and reference images p(r) (top row), weighting function γ(r) (second row), reconstructions using SH 1 regularisation (third row) and STV regularisation (fourth row).
100 % higher compared to the saline. We remark that the reconstructions with the experimental data are only qualitative in the sense that they are computed using 2D models from data that is acquired from a translationally symmetric 3D target. The measurements were carried out with the KIT 4 measurement system [31]. Sixteen adjacent current patterns were used and the voltages differences between the adjacent electrodes were used as measurements, leading to 256 voltage measurements. The amplitude of the injected currents was 5 mA with frequency 10 kHz.
For the estimation of measurement error statistics, 40 000 realizations of the data were measured. The distribution of the measurement noise conformed well with a zero mean Gaussian model. The covariances Γ e,k were formed as sample averages for each of the k = 1, . . . , 16 current patterns separately. In KIT 4, different current patterns use partially different circuit boards and have different switch states, which may result in different measurement error covariances for different current σ true p(r) (partial) p(r) (partial) p(r) (partial) γ(r) SH1 STV Figure 5. Numerical experiment (Case 3): Reconstructions with increasing uncertainty about the edge location in the partial edge information. Rows top to bottom: σ true and reference images p(r) (top row), weighting function γ(r) (second row), reconstructions using SH 1 regularisation (third row) and STV regularisation (fourth row).
patterns. Significant external low frequency sources were not present, and the serial autocorrelation was verified to essentially vanish after the zeroth lag. Thus, the errors in the demodulated voltages could be modelled as mutually independent between current patterns and the overall covariance Γ e was constructed as a block diagonal matrix, Γ e = BlockDiag(Γ e,1 , . . . , Γ e,16 ) The standard deviations (diag(Γ e )) 1/2 of the noise were in the range from 0.01 % to 2.64 % of the mean of the measured voltages and the signal to noise ratio (SNR) was 65.52 dB. The structure of (essentially) non-zero elements in the estimates Γ e,k was tridiagonal, implying that the voltage difference measurements between adjacent electrode pairs that involve electrode are correlated.
The reconstructions were computed using a 2D FEM model. The chest shaped 2D domain Ω was divided to M e = 29 000 triangular elements with M n = 15 233 nodes for approximation of the potential u and to N = 6 187 square pixels for the representation of the conductivity, leading to unknown σ ∈ R 6 187 in the inverse problem.
The electrode contact impedances z = (z 1 , . . . , z 16 ) T ∈ R 16 were estimated before the actual reconstruction by solving a least-squares problem (30) (σ 0 , z) = arg min σ0>0,z>0 where σ 0 ∈ R is a homogeneous conductivity value. The non-linear least-squares problem was minimized using the Gauss-Newton method equipped with a line search algorithm. The line search was implemented using bounded minimization such that positivity conditions are enforced. Similar least-squares approaches for estimation of contact impedances have been proposed in [24,30]. The reconstructed conductivities are shown in figure 7. The top row shows a photograph of the experiment and the reference images p(r). The reference image for the correct edge prior in the third column was obtained by finding boundaries of the agar targets from the photograph of the experiment. The units of the reference image are arbitrary. The reference image for the partial prior in the fourth column was obtained by blurring the reference image of the correct prior using an averaging disk filter. The second row shows the weighting functions γ(r), which in the correct prior case was obtained by (27) and in the partial prior case by (29). The third row shows the reconstructions using the SH 1 regularisation and the fourth row the reconstructions using the STV regularisation.
The results with the experimental data conform qualitatively to the findings of the numerical simulations-the use of structured regularisation leads to significant improvement over the conventional (no structure) regularisation, even in the case of the partial edge information which has reasonably wide uncertainty about the location of the edges.
We remark that the simulated test cases and the experimental test case used different current injections. However, the findings between the different methods were qualitatively similar with both current injections, exemplifying that the proposed approach is not tied to a particular choice of the measurement protocol.
5. Discussion & conclusions. In this paper, we studied incorporation of structural prior information into the electrical impedance tomography problem from an auxiliary reference image using the parallel level sets approach. The approach is based on (general) variational regularisation and it uses a symmetric tensor field for construction of anisotropic weighting for the gradient of the unknown conductivity, given the auxiliary reference image. In this study, we implemented numerically the structured regularisation for quadratic smoothness regularisation and sparsity promoting (smoothed) TV regularisation, and computed reconstructions by minimizing regularised weighted least-squares functional using a lagged Gauss-Newton method. The approach was tested using simulated data and experimental EIT data from a laboratory setup. The results were computed using correct edge information in the prior and partial edge information and compared to the conventional (no structure) regularisation. In the cases studied, the structured regularisation improved the reconstruction accuracy compared to conventional regularisation. Also, in all the cases the structured TV regularisation gave smaller reconstruction errors than the structured smoothness regularisation. While the differences were quite small in the case of correct edge information, they were more pronounced in cases of partial edge information. Further, the structured regularisation was found to improve over the conventional regularisation for both regularisation models even in cases of reasonable uncertainty about the locations of the edges in the prior-this finding is compelling in medical applications, where the structural information could be deliberately blurred to account for organ movement and body deformations, or the model could be constructed based on statistical atlases in the first place. One interesting possibility for constructing the prior in lung monitoring applications could be to use cardiac gated CT and ECG data for estimating the extent of uncertainty in the boundary locations.
In this study the approach was tested using 2D models. The extension to 3D, however, is straightforward. Here the sparsity promoting TV regularisation was based on the smoothed TV functional and the estimates were computed using gradient based optimisation. An extension to non-smooth optimization algorithms with non-linear forward operator is non-trivial and will be a subject of future work.