Colour level set regularization for the electromagnetic imaging of highly discontinuous parameters in 3D

Abstract

In this paper, we propose a novel reconstruction scheme for the low-frequency near-field electromagnetic imaging of high-contrast conductivity distributions inside shielded regions using the system of Maxwell's equations in 3D. In our novel scheme, we focus on estimating the shape characteristics of the electrical conductivity profile inside these regions from low-frequency electromagnetic data measured at external locations for a single frequency. We introduce a colour level set regularization scheme which is a shape-based approach focusing on the simultaneous reconstruction of several shape-like distributions of different conductivity values in the same region of interest. Using two numerical experiments addressing a three-value reconstruction problem related to the imaging of shielded boxes or cargo containers, we compare this novel approach with results obtained from standard voxel-based reconstruction schemes on the one hand and the more established two-value shape-based approach on the other hand. We demonstrate that, depending on the particular situation of the imaging setup, this three-value (or in general multiple-value) shape-based reconstruction technique has the potential to provide superior reconstruction results in many situations, in particular regarding reconstruction of the correct shapes. We also discuss particular challenges of this novel methodology.

Keywords:

• Inverse problems
• electrical conductivity
• shape-based reconstruction
• topological perturbations
• colour level sets

2010 Mathematics Subject Classifications:

• 35R30
• 49N45
• 65F22

1. Introduction

Level set regularization schemes have become a popular methodology when interested in solving shape-based inverse problems. For example, single level set inversion schemes have been designed in a multitude of applications, including non-destructive testing [1–4], reservoir imaging [5–7], medical imaging [8–10], or geophysical prospecting [11,12], amongst others. Its extension, colour level set inversion schemes, have also proven to be useful as they sometimes depict more realistic scenarios. For example, a structural level set method was developed in [13] to detect tumours in breast tissue, a related approach was applied to industrial process tomography in [14], and similar schemes have been developed in the history matching of petroleum reservoirs [15]. Additional level set regularization techniques, such as parametric level sets, have also been used in other imaging modalities such as electrical impedance tomography, see [16–22] for more information.

We focus on developing a scheme to image the contents of boxes and/or containers, using near-field electromagnetic (EM) data. This is a challenging task as the underlying inverse problem is highly ill-posed and the data available is typically limited, resulting in a severely under-determined system. Moreover, shielded containers stop high-frequency EM waves from passing into the contents, through the skin effect. On the other hand, low-frequency contributions to the time signal do have the potential to extract information from inside these shielded regions. The use of low frequency in near-field EM imaging presents an additional challenge of obtaining high resolution reconstructions of material properties in the contents of a shielded container.

In this paper, we introduce a novel colour level set inversion scheme for low-frequency near-field EM imaging to reconstruct material properties inside a shielded region of interest. A special focus is on high-contrast situations where traditional pixel- or voxel-based approaches do not yield satisfactory results due to the inherent smoothing properties of those techniques when using low-frequency data. Those techniques tend to deliver severely low-pass filtered versions of the box content from which it can be impossible to extract sufficient information for detecting or classifying specific objects. In a previous paper [4], a shape-based method as well as a sparsity regularized technique have been proposed in order to better address such a challenging situation. However, whereas both of them perform very well in situations where the content of the domain of interest is mainly composed of two classes of materials with only low variation inside each class, they both struggle to handle other scenarios. For example, when the materials have high contrast with values centred about more than two clearly defined average values. One such situation is the screening of a box or cargo container where there are regions filled with air, as well as moderate and highly conductive material. Classifying such a content into only two classes as implicitly or explicitly done in the sparsity approach or the two-value level set approach does not deliver good results. Therefore, in this paper, we extend the two-value shape reconstruction approach to a multiple-value shape reconstruction approach using more than one level set function for the practical modelling. In particular, we will concentrate on the most prevalent situation of three different classes as indicated in the example given above. This approach is novel in this challenging application and shows surprising results.

The presented colour level set scheme, or its variants, are not limited to low-frequency electromagnetic applications but can be applied to a wide range of tomography problems. This includes X-ray computerized tomography (CT), electrical impedance tomography, microwave medical imaging, and history matching of production data in reservoir characterization. Some currently available work on several of those applications can be found in [13,15,17,20,23,24], and in the references provided there.

It is interesting to compare our application of low-frequency electromagnetic imaging with the above mentioned application of CT, which essentially uses very high-frequency electromagnetic beams (X-rays). In most clinical applications of CT (nowadays using sophisticated sampling geometries) the classical reconstruction algorithms provide high-resolution images [25], such that the combined reconstruction and segmentation approach as followed in our discussion does not on the first view seem to provide significant advantages there. Those CT reconstructions are usually obtained by applying some approximation of the Radon inversion formula, which result in highly efficient implementations, provided that the raw data (so-called ‘sinograms’ or photon counts) are sufficiently complete and of sufficiently high quality (see [25] for details). Therefore, when segmented CT images are needed, a standard CT image is obtained using fast reconstruction techniques, which is then further treated by tailor-made (for example, level set based) segmentation post-processing techniques in order to obtain segmented images [26].

However, even in X-ray CT there are many situations where those quality conditions on raw data are not satisfied. This is the case for example in many industrial settings, but also in specific clinical and biochemical applications where only limited view data or very sparse data are available or desired. In those situations, a combined level set reconstruction and segmentation algorithm, similar to the one proposed here, becomes much more competitive or even superior to standard reconstruction techniques since many classical reconstruction algorithms show significant artefacts which cannot easily be handled by standard segmentation post-processing strategies. For some discussions in the literature addressing such situations we refer to [10,25,27–30].

In the application of low-frequency electromagnetic imaging, there currently does not exist the luxury of a high-resolution reconstruction algorithm (as in X-ray CT) in the first place. Compared to standard X-ray CT images, reconstructions show significantly less details here and interfaces which might be present in the original objects are usually not at all represented in the images obtained from standard reconstruction techniques [31,32]. They need to be enforced by tailor-made reconstruction techniques, making sure that they are in agreement with the physically measured data. The use of a combined ‘reconstruction and segmentation’ algorithm as discussed in this paper seems to be a very natural way of obtaining such segmented images from physical data.

The remaining part of the paper is organized as follows: Section 2 describes the mathematical setup of our low-frequency EM imaging problem with Maxwell's equations. Section 3 formulates the underlying EM inverse problem for recovery of the conductivity value inside the imaging domain and derives some general concepts which are essential for voxel-based and colour level set based reconstruction schemes. As a side product, a voxel-based reconstruction scheme is derived in this section which we use for comparison with the colour level set scheme. Section 4 introduces the colour level set shape-based scheme and formulates the corresponding inverse problem. Section 5 presents a methodology for avoiding certain local minima in the colour level set reconstruction. Section 6 presents two numerical experiments for the colour level set inversion as well as a comparison with alternative regularization schemes (including results using the methodology in Section 5). In this section, we also provide reconstructions which are based on more realistic compositions inside the region of interest in the true phantom with significant internal parameter fluctuations in one of the target zones. Section 7 presents a summary and conclusions of our findings. Lastly, the Appendix explains in details the new tailor-made line search method for our scheme.

2. The near field electromagnetic imaging problem

We model the propagation of EM fields by Maxwell's equations in the frequency domain. These are given as: (1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) (1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) where we assume a time-dependence $\exp (-i\omega t)$ for a given and fixed frequency $\omega =2\pi f$, with $i=\sqrt{-1}$ denoting the imaginary unit and t denoting time. In (1a(1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) ), (1b(1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) ), $\mathbf{E}$ denotes the electric field, $\mathbf{H}$ the magnetic field, $\mathbf{M}$ the magnetic source and $\mathbf{J}$ the electric source. The material parameters are (2) $a(\mathbf{x})=i\omega \mu (\mathbf{x}),b(\mathbf{x})=\sigma (\mathbf{x})-i\omega ϵ(\mathbf{x}),$(2) where $\sigma (\mathbf{x})$ is the electrical conductivity, $\mu (\mathbf{x})$ is the magnetic permeability and $ϵ(\mathbf{x})$ is the electric permittivity. The index j in (1a(1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) ), (1b(1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) ) indicates that the jth source distribution is considered, where $j=0,1,\dots ,n_s-1$. Therefore, in our setup, we will consider $n_s$ different applied sources. These will be modelled in our numerical experiments as rectangular wire loops, but the reconstruction schemes presented here do not depend on this particular choice. The sources give rise to probing fields ${\mathbf{E}}_j$ and ${\mathbf{H}}_j$ in (1a(1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) ), (1b(1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) ), which can be measured at the receiver locations. Again, we will use rectangular wire loops as receiving devices, even though the general schemes derived here will work with other measurement approaches with only very minor adjustments. The physically obtained (or in our numerical experiments simulated) measurements are called in the following simply ‘data’, which will be used in order to infer medium characteristics at those places which are not accessible by the antennas.

We are mainly interested here in the non-invasive imaging of the interior of boxes whose walls show a significant conductivity profile. The main difficulty is that EM fields of high frequency do not penetrate sufficiently through such conductive walls due to the well-known skin-effect. Therefore, relatively low frequencies need to be employed which require specialized approaches for data inversion. A possible experimental setup, which is used in our numerical experiments, is indicated in Figure 2. It also shows a possible distribution of sources and receivers outside the shielding walls. Here, the domain of interest is represented by a cuboid-shaped domain $\mathrm{\Omega }\subset {\mathbb{R}}^3$ surrounded at all six sides by shielding walls of constant and known thickness. The sources and receivers are distributed uniformly along two planes parallel to the x−y plane each of constant z coordinate. This means that they are opposing each other, separated by the shielded box of interest (see Figure 2). Notice that there are no sources or receivers at the other four sides of the box, which limits the aperture of our measurement setup but is more realistic in potential applications such as the fast screening of luggage or cargo containers. We mention again that most of these special assumptions on the geometry can be relaxed in order to arrive at more general imaging scenarios but are sufficient for our proof-of-concept study presented here.

In this paper, we are interested in estimating shapes of the electrical conductivity profile $\sigma (\mathbf{x})=\mathrm{\Re }b(\mathbf{x})$ (where ℜ refers to taking the real part of the complex quantity) inside the box from data measured outside the shields. For simplicity, we choose $\mu (\mathbf{x})={\mu }_0$ and $ϵ(\mathbf{x})={ϵ}_0$ everywhere where ${\mu }_0$ and ${ϵ}_0$ represent free space values. This means two of the three material parameters are constant and known throughout the entire domain. We will therefore often refer to the parameter b of (2(2) $a(\mathbf{x})=i\omega \mu (\mathbf{x}),b(\mathbf{x})=\sigma (\mathbf{x})-i\omega ϵ(\mathbf{x}),$(2) ) simply as ‘conductivity’ in this paper, even though it also contains the known parameter ${ϵ}_0$. In our future research, we plan to address the incorporation of two or three material parameters as space-varying unknowns into the inversion process, but in this study, we will concentrate on the conductivity σ only.

We require a mathematical model of the scenery when doing the inversion. Traditionally, a discretization of the entire 3D domain into a (often rectangular) grid of voxels is the most popular approach for the inversion, where each voxel is assigned a potentially different value of conductivity inferred from data in the inversion process. Following notation from the classical (mainly 2D) literature, we will call this approach ‘pixel-based inversion’, even though it is actually voxels that are assigned values to in our 3D application. In our approach, however, we have chosen to use a completely different model for reasons explained further below. We are interested in imaging situations where some objects of significantly higher conductivity profile are located inside the box surrounded by some background of lower conductivity values, and the task is to identify and characterize those high-conductivity objects from data obtained outside the box. In pixel-based inversion approaches, the classification into different material types is done as a post-processing step, after an inversion has been obtained, by some form of segmentation routine. However, a potential problem occurs here since the classical segmentation post-processing tools, being purely image based but not data based, do not usually verify the agreement of the segmented image with the physically obtained data. Since the modifications which are needed by classical segmentation tools for transforming a very smooth reconstruction into a segmented image are significant in our situations, there is no guarantee that the segmented image still agrees with the measured data the same way as the originally smooth reconstruction.

Therefore, in a previous publication [4], we have proposed a single-level set reconstruction approach which combines the segmentation task and the inversion task into one optimization loop, providing segmented 3D images honouring the physically measured data. In that approach, it is assumed that only two conductivity values are permitted inside the box, namely the high-conductivity value of the object of interest and one single lower value for the background medium. This approach yields reasonable results in many applications. However, there are certain situations where the background medium itself is highly heterogeneous (e.g. consisting of air and some lower conductive materials), such that modelling it by a constant average value everywhere seems to be an oversimplification. In this paper, we extend the previous model to also allow the background medium to obtain two (or even more) different (fixed and suitably chosen) values of a priori unknown distribution. This means, that now in total three (or a specified low number of) different conductivity values are considered for each location inside the box, one for the object of interest and the remaining values for the background. In the inverse problem, we are trying to identify shapes of not too complicated topology which represent distributions of the different materials in the box. During the inversion, each background conductivity value will be classified this way, by implicit means, to belong to exactly one of the corresponding prescribed background material classes. Embedded in this segmented background is the object of interest which corresponds to a separate class representing high-conductivity values. In other words, we perform a simultaneous reconstruction and segmentation of the domain of interest into various different zones which, when we plug the corresponding distribution of material parameters into the forward simulator, reproduces the measured data.

We mention again that this classification can be derived and extended in a flexible way to represent an arbitrary (usually small) number of representative classes. See for example a previous publication related to a different application [13] where a medical setup of higher (microwave) frequency ranges has been treated successfully with more than three parameter classes using a similar scheme, but only in 2D so far. We have tested the use of more than three classes in our numerical experiments (not presented in this paper) also to the low-frequency situation discussed in this paper but have observed that the use of more than three classes does not provide any advantages at this low-frequency regime, mainly due to the lack of resolution inherent to it. Therefore, we will concentrate in this paper on the particular setup of three material classes only. It will become obvious below how to extend our model to incorporate more than three classes in a straightforward way if desired.

In our numerical experiments, we assign a very low value $\sigma \approx {10}^{-8}$ Sm${}^{-1}$ to the region which represents poorly conductive materials such as air or other essentially non-conductive materials. In addition to this very low conductivity region, there will also be one moderate value representing a class of material at the intermediate regime, and one highly conductive region capturing any materials with significant conductivity properties, all of unknown shape. We emphasize that in our experience it is not very critical which values exactly are taken to represent those classes due to the significant contrast between them. However, in cases where prior information is available indicating the use of specific values this can and should certainly be incorporated in the algorithm by assigning those numbers to the individual classes. A summary of conductivities used in our numerical experiments is described in Table 1.

Table 1. Conductivity parameters used in numerical experiments.

Domain characteristic	$\sigma =\mathrm{\Re }(b)$ (Sm${}^{-1}$)
Region $S_3$	$1\times {10}^{-8}$
Region $S_2$	1.0 / 10.0
Region $S_1$	0.5
Shields	0.1

The shielded walls, which are not part of our level set representation and are modelled explicitly, are here assumed to have a known topology and a moderate conductivity value of $\sigma \approx 0.1$ Sm${}^{-1}$. This assumption will be relaxed in our future work when we incorporate estimation of some characteristics of the shields into the inversion task. An important observation of Table 1 is that the conductivity values range over several orders of magnitude. This is difficult to model with a pure pixel-based inversion model, due to the typical effect of smeared-out values in the reconstructions. In those classical approaches, sharp interfaces between different materials are usually smeared out over large transition zones, which makes the identification of particular regions or shapes extremely difficult or impossible. We will demonstrate this in our numerical experiments.

3. Pixel-based inversion tools

Inverse problems for Maxwell's equations have been investigated for a long time in a variety of applications and settings. For recent overviews of available theoretical and computational results, we refer to [33–35] and the references provided there. The important questions of uniqueness and differentiability of related inverse problems are addressed for example in [36–41]. For derivation of the algorithms in this paper, we assume differentiability (in properly chosen function spaces) of all forward mappings that arise here, allowing us to adopt well-established expressions for the relevant derivatives (or sensitivities) with respect to the unknown medium parameters [35,39,42–46].

The colour level set inversion approach outlined above can be considered a shape-based regularization scheme which provides an alternative to the more classical Tikhonov-regularized pixel-based inversion schemes. As already mentioned, due to the smearing effect related to Tikhonov regularization, the classical Tikhonov-regularization pixel-based scheme is not well-suited for the imaging of high-contrast situations with low frequency EM fields. Nevertheless, both approaches are based on gradient calculations with respect to the continuous parameter-to-data maps, which need to be considered first. Fortunately, most of the results related to this basic ingredient are already well-known, see for example [43,47] or some of the above cited overview articles. Therefore, we will very briefly summarize some of the most important results. These gradient calculations are relevant to both pixel-based and shape-based approaches. They will also provide us directly with a basic Kaczmarz-type pixel-based inversion approach which we will use in our numerical experiments for comparing our novel shape-based techniques with more classical Tikhonov-style schemes.

Let us index the forward operator by each source j, giving rise to the notation ${\mathcal{F}}_j[b]={\mathcal{M}}_j{\mathbf{E}}_j(\mathbf{x};b)$. The quantities ${\mathbf{E}}_j$, ${\mathbf{H}}_j$ solve (1a(1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) ), (1b(1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) ) with parameter b and ${\mathcal{M}}_j$ is a linear measurement operator which might be dependent on the source. Here we take ${\mathcal{M}}_j$ to behave the same across all sources. In our numerical setup explained further below and in [48], we will formally model measurements at each of the $n_r$ rectangular closed-loop receiver coils $S_k$ ($k=0,\dots ,n_r-1$) by (3) ${\left[{\mathcal{M}}_j{\mathbf{E}}_j\right]}_k={\oint }_{S_k}{\mathbf{E}}_j\cdot \hat{\mathbf{t}}\mathrm{d}l,$(3) assuming proper regularity of ${\mathbf{E}}_j$. Here, $\hat{\mathbf{t}}$ denotes the tangential unit vector along such a receiver coil. More general linear measurement operators can be used as well without significant modifications of the algorithm. See for example [32,39,44,48,49] for more details and for alternative formulations.

A typical way of addressing the pixel-based inverse problem is to write down a least squares data misfit functional in terms of residual operators. In our case this reads as (4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) with ${\mathcal{R}}_j[b]={\mathcal{F}}_j[b]-{\tilde{\mathcal{F}}}_j$, and $\Vert \cdot {\Vert }_{Z_j}$ being the norm in the data space $Z_j$. $\widetilde{{\mathcal{F}}_j}$ indicates the true data. Here, the index j refers again to the chosen source distribution. Therefore, the residual operator quantifies the mismatch between model prediction (forward problem) and true data. The cost functional $\mathcal{J}[b]$ is useful for monitoring progress of any iterative estimation technique and is often used directly for the design of suitable reconstruction schemes. For example, we can formulate an optimization problem incorporating $\mathcal{J}[b]$ with or without additional suitably chosen regularization terms. Many popular minimization schemes use the formal gradient of (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) ): (5) $\mathrm{\nabla }\mathcal{J}[b]=\sum _{j=0}^{n_s-1}\mathrm{\nabla }{\mathcal{J}}_j[b],$(5) which can be calculated efficiently by a so-called adjoint scheme [42,43]. We refer the reader to these references for a full derivation of the adjoint problem and the gradient of (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) ). Here, it suffices to quote the result. This is found to be: (6) $\mathrm{\nabla }{\mathcal{J}}_j[b]={\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]$(6) and is given by (7) $[{\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]](\mathbf{x})=\overline{{\mathbf{E}}_j(\mathbf{x})\cdot {\mathcal{E}}_j(\mathbf{x})}.$(7) Here, the operator $R_j^{\prime }[b{]}^{\ast }$ is known as the adjoint linearized residual operator and the overline on the right hand side indicates to take the complex conjugate. The fields ${\mathbf{E}}_j(\mathbf{x})$ are taken from (1a(1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) ), (1b(1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) ) and ${\mathcal{E}}_j(\mathbf{x})$ are obtained by solving a supplementary adjoint Maxwell problem of similar structure as (1a(1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) ), (1b(1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) ) but with using the residual values ${\mathcal{R}}_j[b]$ as artificial adjoint sources at the receiver locations. For more information on this scheme, see [4,35,42,43,45,46,48].

The non-linear problem described is of large scale in 3D. Usually, iterative techniques are employed for the solution, requiring repeated calculation of descent directions to $\mathcal{J}[b]$ in (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) ). Many standard schemes require calculation of the full gradient $\mathrm{\nabla }\mathcal{J}[b]$ in (5(5) $\mathrm{\nabla }\mathcal{J}[b]=\sum _{j=0}^{n_s-1}\mathrm{\nabla }{\mathcal{J}}_j[b],$(5) ). Depending on the available Maxwell solver, this might consume considerable resources and processing time in each iteration. Therefore, alternative techniques have been developed in order to speed up the inversion process. One of those is the non-linear Landweber-Kaczmarz (LK) scheme which is employed here. The optimizer cycles over each source in a certain order and only requires to calculate gradients of partial data sets in each step, which can be done efficiently with the above mentioned adjoint method.

With this single-step scheme it is no longer the primary goal to maximally reduce the data misfit $\mathcal{J}[b]$ in (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) ) in each individual step. Instead, in each iteration, we reduce the data misfit ${\mathcal{J}}_j[b]$ with respect to the jth source only, which needs to be done very carefully in order not to deteriorate too much the data fit regarding the remaining sources. For more information on this general scheme, we refer to [28,42,50–53].

The scheme yields a single-step update formula: (8) $b^{s+1}=b^s-\tau \mathrm{\Re }(\mathrm{\nabla }{\mathcal{J}}_{[s]}[b^s]),$(8) where we will follow here a sequential rule $[s]:=(s\text{modulo}{n}_s)\in \{0,1,\dots ,n_s-1\}$ for simplicity. Notice that we take the real part ℜ in (8(8) $b^{s+1}=b^s-\tau \mathrm{\Re }(\mathrm{\nabla }{\mathcal{J}}_{[s]}[b^s]),$(8) ) since we assume that the electric permittivity is known. As already mentioned, choosing the step size τ in this single-step scheme is a difficult problem, and it should be chosen such that updates maintain a balance with the already achieved data fit regarding the remaining sources. Convergence properties of such a single-step approach have been discussed in the literature, see for example [28,52]. A random selection of source positions, leading to variants of stochastic gradient schemes, might provide additional benefits, see for example [54] for related applications in the field of machine learning.

The iteration formula in (8(8) $b^{s+1}=b^s-\tau \mathrm{\Re }(\mathrm{\nabla }{\mathcal{J}}_{[s]}[b^s]),$(8) ) updates the electrical conductivity profile such that it minimizes the cost functional given in (5(5) $\mathrm{\nabla }\mathcal{J}[b]=\sum _{j=0}^{n_s-1}\mathrm{\nabla }{\mathcal{J}}_j[b],$(5) ). This minimization usually delivers quite smooth results due to the specific form of individual updates in (8(8) $b^{s+1}=b^s-\tau \mathrm{\Re }(\mathrm{\nabla }{\mathcal{J}}_{[s]}[b^s]),$(8) ) but also can contain quite erratic features since in its current form, there are no additional regularization terms, implicit or explicit, included. Therefore, classical regularization approaches often add a Tikhonov-type additional term to (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) ) in order to stabilize the inversion process. This is usually called ‘Tikhonov regularization’. A similar effect is achieved by keeping (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) ) unchanged but instead using specifically chosen function spaces. In particular the choice of certain Sobolev spaces enforces higher regularity properties of their members. We will adopt such a Sobolev gradient approach because it does not come with the need for modifying the cost functional (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) ). In particular, a direct comparison is possible between the results obtained with this approach and those obtained by our colour level set shape reconstruction method, since both are designed to minimize the same data misfit functional but with respect to a different set of optimization parameters.

In the Sobolev gradient approach, we assume that the quantities of interest (say, u, v) are members of a Sobolev space $W_{\alpha ,\beta }$ equipped with the $(\alpha ,\beta )$ norm $\Vert u{\Vert }_{\alpha ,\beta }^2=⟨u,u{⟩}_{\alpha ,\beta }$ and inner product $⟨u,v{⟩}_{\alpha ,\beta }=\alpha ⟨u,v{⟩}_{L^2}+\beta ⟨\mathrm{\nabla }u,\mathrm{\nabla }v{⟩}_{L^2}$ with suitable boundary conditions, see for example [3,4,9,11,24] for more details. In our numerical experiments, we will choose $\alpha =1$ and select $0<\beta <1$ empirically in order to achieve a stable reconstruction process. Since the gradient operator depends on the spaces and corresponding inner products, we have to calculate new expressions for the corresponding gradient directions. Fortunately, as shown in [3,4,9,11,24], the Sobolev gradient expressions ${\widehat{\mathrm{\nabla }\mathcal{J}}}_j[b]$ can be obtained in a straighforward way from the ${\mathcal{L}}^2$ based gradients (6(6) $\mathrm{\nabla }{\mathcal{J}}_j[b]={\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]$(6) ), (7(7) $[{\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]](\mathbf{x})=\overline{{\mathbf{E}}_j(\mathbf{x})\cdot {\mathcal{E}}_j(\mathbf{x})}.$(7) ) by simply projecting these expressions into a smoother space using the following recipe (9) $(\alpha I-\beta \mathrm{\Delta }){\widehat{\mathrm{\nabla }\mathcal{J}}}_j[b]=\mathrm{\nabla }{\mathcal{J}}_j[b]\text{in}\mathrm{\Omega }.$(9) Here, I is the identity operator and Δ is the Laplacian operator. Solving (9(9) $(\alpha I-\beta \mathrm{\Delta }){\widehat{\mathrm{\nabla }\mathcal{J}}}_j[b]=\mathrm{\nabla }{\mathcal{J}}_j[b]\text{in}\mathrm{\Omega }.$(9) ) to obtain ${\widehat{\mathrm{\nabla }\mathcal{J}}}_j[b]$ provides us with a smoothed gradient for the inverse problem. For derivation of (9(9) $(\alpha I-\beta \mathrm{\Delta }){\widehat{\mathrm{\nabla }\mathcal{J}}}_j[b]=\mathrm{\nabla }{\mathcal{J}}_j[b]\text{in}\mathrm{\Omega }.$(9) ) and for more information on practical ways of calculating this smoothed gradient, we refer the reader to a more detailed discussion in [3,4,9,11,24].

4. Colour level set shape-based inversion

Level set shape-based reconstruction methods provide an alternative to pixel-based schemes and have been applied in many areas [23,24,47,55–57]. Level set methods provide additional regularization to the inverse problem as they combine parameter values into a small number of classes with fixed (or simultaneously estimated) representative values. This way level set based schemes solve two tasks at the same time: reconstruction of physical profiles from physically measured indirect data and the segmentation of that profile into a small number of classes. The most basic of those schemes, here called standard level set method, assumes that the domain is composed of exactly two different classes of fixed parameter values. They have been shown to perform well in certain situations of low-frequency EM inverse problems [1,4,7,11,23,24]. However, in some scenarios, for example, when multiple conductivity values of significantly different magnitudes are present in the domain, standard level set methods are not optimal (neither are pixel-based inversion schemes). In such scenarios, colour level set methods can be applied in an attempt to deal with such more challenging situations.

We refer to the theory of standard level set methods, in particular with respect to the application considered here, to a previous publication [4]. In the following, we directly turn our attention on the extended setup necessary for handling more than two different regions by a certain number of level set functions. A very general introduction to multiple shape reconstruction by a variety of practical level set based approaches can also be found in [23,24].

For modelling three different regions, we will be using here two different level set functions. Therefore, we begin by introducing two sufficiently smooth level set functions ${\varphi }_{1,2}:{\mathbb{R}}^3\to \mathbb{R}$, such that (10) $b(\mathrm{\Phi })(\mathbf{x})=\{\begin{array}{ll}{b}^1& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})>0,\end{array}$(10) where $\mathrm{\Phi }=({\varphi }_1,{\varphi }_2)$ and $(b^1,b^2,b^3)$ are the corresponding parameter values (the superscripts not to be confused in the following with powers of b which are not needed in this study).

A similar but conceptually different multiple level set model has been introduced originally for the application of image segmentation in [58] where two level set functions are used for describing four different regions. In contrast to the model used in [58], our representation (10(10) $b(\mathrm{\Phi })(\mathbf{x})=\{\begin{array}{ll}{b}^1& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})>0,\end{array}$(10) ) uses two level set functions in order to only model three different regions. Moreover, our model is not completely symmetric in the three values $(b^1,b^2,b^3)$ in the sense that ${\varphi }_1$ exclusively specifies the region filled with $b^1$, whereas the sign of ${\varphi }_2$ by itself does not provide sufficient information about which value to choose. See Figure 3 for a visualization of this situation. This hierarchical arrangement actually provides some advantages compared to the model proposed in [58] when applied to shape reconstruction problems from indirectly obtained data.

One of its advantages is that it allows us in a very convenient way to introduce single new regions at any stage of the reconstruction process into the model without the need of redefining the representation of already existing regions. We just add one additional level set function which will specify a new region without the risk of ambiguity. For example, a reconstruction could start from just using ${\varphi }_1$ during early sweeps, segmenting the domain into only two regions, and at a later stage could branch out by adding ${\varphi }_2$ in order to refine the segmentation.

Another advantage is that two regions usually do not have a point-like (or cusp-like) interface with each other (except of some very unlikely special situations). Instead, two different regions normally share a full curve-like interface in 2D or a full surface-like interface in 3D, which makes it easier to apply the concept of shape velocity and shape gradient for shape propagation or their ‘narrowband analogon’ as it is outlined further below and used in this paper. Notice also that our scheme extends accordingly when more than two level set functions are used in order to model more than three regions and avoids some local minima which can occur with some other schemes of using multiple level set functions [23,24].

In general, all three profiles $(b^1,b^2,b^3)$ can be either smoothly varying functions in these regions or can be taken as constants. In our numerical experiments and in the notation of this paper, we only consider the latter choice. An equivalent way of writing the conductivity profile in (10(10) $b(\mathrm{\Phi })(\mathbf{x})=\{\begin{array}{ll}{b}^1& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})>0,\end{array}$(10) ) is (11) $b(\mathrm{\Phi })=b^1(1-H({\varphi }_1))+b^{2}H({\varphi }_1)[1-H({\varphi }_2)]+b^{3}H({\varphi }_1)H({\varphi }_2),$(11) where $H(\cdot )$ is the Heaviside function.

With the colour level set representation of the electrical conductivity profile formally defined, we now introduce the relevant data misfit and cost functional: (12) $\mathcal{T}(\mathrm{\Phi })=\mathcal{R}[b(\mathrm{\Phi })],\mathcal{J}(\mathrm{\Phi })=||\mathcal{T}(\mathrm{\Phi })|{|}_{{\mathcal{L}}^2(\mathrm{\Omega })}^2,$(12) (and similarly the quantities ${\mathcal{T}}_j(\mathrm{\Phi })$ and ${\mathcal{J}}_j(\mathrm{\Phi })$ as in (4(4) $\mathcal{J}[b]=\sum _{j=0}^{n_s-1}{\mathcal{J}}_j[b];\text{and}{\mathcal{J}}_j[b]:=\frac{1}{2}{\left||{\mathcal{R}}_j[b]|\right|}_{Z_j}^2,$(4) )). We are now interested in finding a suitable $\mathbf{\Phi }$ such that the data misfit $\mathcal{J}(\mathrm{\Phi })$ is minimized. Notice that, even though using the same symbol as in the pixel-based case for simplicity in the notation, $\mathcal{J}$ is now assumed to be a function of $\mathbf{\Phi }$.

When comparing with more traditional colour level set approaches for image processing, we would like to point out here that the cost functional (12(12) $\mathcal{T}(\mathrm{\Phi })=\mathcal{R}[b(\mathrm{\Phi })],\mathcal{J}(\mathrm{\Phi })=||\mathcal{T}(\mathrm{\Phi })|{|}_{{\mathcal{L}}^2(\mathrm{\Omega })}^2,$(12) ) for our colour level set inversion method still contains the data-to-model mapping $\mathcal{R}$ in its formulation. Image post-processing applications usually replace this (physical) data fidelity term by a term measuring proximity to a reference image [26,58] as outlined in Figure 1. We cannot rely on any reference image in our application but only on the physically gathered data.

Figure 1. Comparison of the traditional 'first reconstruct then segment' scheme with the proposed 'simultaneously reconstruct and segment' scheme. The results on the right hand side are both segmented images, of which only the one from the simultaneous approach is guaranteed to satisfy the physically obtained EM data.

With (5(5) $\mathrm{\nabla }\mathcal{J}[b]=\sum _{j=0}^{n_s-1}\mathrm{\nabla }{\mathcal{J}}_j[b],$(5) ) and (6(6) $\mathrm{\nabla }{\mathcal{J}}_j[b]={\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]$(6) ), and by applying formally the chain rule for differentiation with respect to $\mathbf{\Phi }$, we obtain an expression for the gradient of the colour level set cost functional: (13) $\mathrm{\nabla }\mathcal{J}(\mathrm{\Phi })=\sum _{j=0}^{n_s-1}\mathrm{\nabla }{\mathcal{J}}_j(\mathrm{\Phi })=\sum _{j=0}^{n_s-1}{\mathcal{T}}_j^{\prime }(\mathrm{\Phi }{)}^{\ast }{\mathcal{T}}_j(\mathrm{\Phi }),$(13) where (14) ${\mathcal{T}}_j^{\prime }(\mathrm{\Phi }{)}^{\ast }{\mathcal{T}}_j(\mathrm{\Phi })=\mathrm{\Re }\left[\begin{array}{c}{b}_{{\varphi }_1}(\mathrm{\Phi }){\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]\\ b_{{\varphi }_2}(\mathrm{\Phi }){\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]\end{array}\right],$(14) and (15) $\begin{array}{rl}{b}_{{\varphi }_1}(\mathrm{\Phi })& =[(b^2-b^1)+(b^3-b^2)H({\varphi }_2)]\delta ({\varphi }_1),\end{array}$(15) (16) $\begin{array}{rl}{b}_{{\varphi }_2}(\mathrm{\Phi })& =(b^3-b^2)H({\varphi }_1)\delta ({\varphi }_2).\end{array}$(16) Here, the notation $b_{{\varphi }_i}$ refers to the partial derivative with respect to that level set function, ℜ indicates again to take the real part in order to account for real-valued level set functions, and $\delta (\cdot )$ denotes the Dirac delta distribution. Evaluating those expressions (14(14) ${\mathcal{T}}_j^{\prime }(\mathrm{\Phi }{)}^{\ast }{\mathcal{T}}_j(\mathrm{\Phi })=\mathrm{\Re }\left[\begin{array}{c}{b}_{{\varphi }_1}(\mathrm{\Phi }){\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]\\ b_{{\varphi }_2}(\mathrm{\Phi }){\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]\end{array}\right],$(14) ) requires us to solve one adjoint Maxwell problem for each index j in order to obtain ${\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]$. In the advent of numerical experiments, $\delta ({\varphi }_i)$ has to be approximated by a function in a suitable space, since it is not well defined on a numerical grid. Therefore, we choose to use already in this theoretical derivation a narrowband function as an approximation to it [59,60], such that $\delta ({\varphi }_i)\approx {\chi }_{B_d({\mathrm{\Gamma }}_i)}$, where (17) $B_{d_i}({\mathrm{\Gamma }}_i)=\{\mathbf{x}:\Vert \mathbf{x}-{\mathrm{\Gamma }}_i{\Vert }_2^2\le d_i/2\},$(17) $d_i>0$ and ${\chi }_D$ is the characteristic function. Notice that this approximation conserves the descent direction property of the (negative) gradient for small values of $d_i$, which is all we need for implementing our descent scheme. In the general case, we assume that boundaries of objects are smooth. Therefore, each component in the gradient from (14(14) ${\mathcal{T}}_j^{\prime }(\mathrm{\Phi }{)}^{\ast }{\mathcal{T}}_j(\mathrm{\Phi })=\mathrm{\Re }\left[\begin{array}{c}{b}_{{\varphi }_1}(\mathrm{\Phi }){\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]\\ b_{{\varphi }_2}(\mathrm{\Phi }){\mathcal{R}}_j^{\prime }[b{]}^{\ast }{\mathcal{R}}_j[b]\end{array}\right],$(14) ) is projected to a smoother space using the recipe in (9(9) $(\alpha I-\beta \mathrm{\Delta }){\widehat{\mathrm{\nabla }\mathcal{J}}}_j[b]=\mathrm{\nabla }{\mathcal{J}}_j[b]\text{in}\mathrm{\Omega }.$(9) ) and is denoted accordingly by ${\widehat{\mathrm{\nabla }\mathcal{J}}}_{[s]}$. The resulting gradient is now a sufficiently smooth function well-defined on the entire domain Ω which can be used for updating the level set function $\mathbf{\Phi }$ everywhere. We mention that the process of extending the Dirac delta to a smooth function on the entire domain is reminiscent of finding an ‘extension velocity’ for classical level set applications in computational physics [60]. Indeed, a similar interpretation can be obtained in our case [23,24]. This leads to the LK-single-step iteration formula: (18) ${\mathrm{\Phi }}_{j+1}^s={\mathrm{\Phi }}_j^s-{\mathit{\tau }}^s\circ {\widehat{\mathrm{\nabla }\mathcal{J}}}_j({\mathrm{\Phi }}_j^s),$(18) where the step size is ${\mathit{\tau }}^s=[{\tau }_1^s,{\tau }_2^s{]}^T$ and $\circ$ denotes the Hadamard (element-wise) product. When using non-linear LK schemes in single-level set inverse problems, choosing the step size has proven to be tricky [3,4]. However the line search scheme devised in [4] was shown to be successful, therefore, we adopt the criterion developed there and extend it to the colour level set regime. Note that the parameter s denotes the current sweep number, which increases by once when we have computed an update for all sources.

The line search essentially uses a backtracking scheme which does not require any additional solves of the forward or adjoint problem. It controls the number of interior voxels whose conductivities change per iteration. A sufficiently large line search parameter τ is chosen and reduced until a smooth evolution of the shape is guaranteed. This amounts to counting the number of voxels that change value based on the current step size, then reducing it until the number of voxels that change value is inside an admissible interval. An outline of the scheme in a generalized colour level set setting is given in Appendix, where we apply the criterion for a particular case.

5. A stochastic seeding process

We mentioned above that the colour level set scheme  (10(10) $b(\mathrm{\Phi })(\mathbf{x})=\{\begin{array}{ll}{b}^1& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})>0,\end{array}$(10) ), (11(11) $b(\mathrm{\Phi })=b^1(1-H({\varphi }_1))+b^{2}H({\varphi }_1)[1-H({\varphi }_2)]+b^{3}H({\varphi }_1)H({\varphi }_2),$(11) ) helps avoiding certain local minima. However, shape evolution schemes with more than two values still suffer from a particular type of local minima due to the nature of only moving interfaces locally.

In particular, certain local minima are associated with the incorrect classification of isolated objects embedded in a background region $S_3$ as being composed of material associated with region $S_1$ instead of region $S_2$, or vice versa. This situation does not occur in standard level set applications where only one background profile and one possible inclusion value are used. We will demonstrate this ambiguity in our numerical experiments.

In order to deal with that type of local minima, we introduce a stochastic seeding process which improves the situation to a certain degree. Unfortunately, due to the high degree of ill-posedness of this inverse problem, and the poor resolution resulting from it, this problem cannot be completely avoided.

The general idea is to randomly place small objects, alternating between type $S_1$ and type $S_2$, into the domain in order to let the algorithm decide whether such an object can favourably replace the already existing structures. Therefore, this strategy can be considered as a way of probing at different stages of the algorithm different variants of ‘advanced initial guesses’, with the hope that the algorithm will select the one leading to the global minimum. Since the entire shape evolution does take a significant number of iterations also without this seeding, such an embedded probing is more efficient compared to starting the entire reconstruction algorithm with a large number of different initial guesses in an ensemble type approach [61]. It still adds computational expenses to the reconstruction, and we will investigate in this paper amongst others whether this additional cost is rewarded by the results or not. For similar stochastic seeding algorithms in a completely different level set reconstruction problem, and employed for different reasons, we refer the reader to the publications [5,6,15,62]. In the following, we very briefly outline the general approach including some technical aspects of it.

Let us assume that we want to place an object centred in a small region ${\mathrm{\Omega }}_2\subset \mathrm{\Omega }$. This amounts to lowering the level set function until it is negative in the desired region, forming a new object. Since our numerical experiments use a much coarser grid compared to [5,6,15], we will only consider two regions: ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$, such that $\mathrm{\Omega }={\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2$, where ${\mathrm{\Omega }}_2$ is the domain containing a seeded object and ${\mathrm{\Omega }}_1$ is its complement. Whilst an additional interim region of smoothing the level set function between ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ would be useful, in our application we will not consider this additional region due to the coarseness of the grid in 3D. Placing an object in $\mathbf{\Omega }$ is practically achieved by solving the following supplementary minimization problem on the corresponding level set function φ at selected steps of the colour level set reconstruction: (19) ${\varphi }_{new}=\underset{\mathrm{\Psi }}{\arg \min }K(\mathrm{\Psi })=\frac{{\alpha }_1}{2}||(\mathrm{\Psi }-{\varphi }_{old}){\chi }_1|{|}^2+\frac{{\alpha }_2}{2}||(\mathrm{\Psi }-\mu ){\chi }_2|{|}^2,$(19) where ${\alpha }_1,{\alpha }_2>0$ are suitably chosen weighting parameters, $\mu <0$ and ${\chi }_2$ and ${\chi }_1$ are characteristic functions concentrated on the domains ${\mathrm{\Omega }}_2$ and ${\mathrm{\Omega }}_1=\mathrm{\Omega }\mathrm{\setminus }{\mathrm{\Omega }}_2$, respectively. The minimizer of this cost functional, ${\varphi }_{new}$, replaces the level set function ${\varphi }_{old}$. The first term in the cost functional is a penalty on the distance between new and old level set functions in the domain outside the seeded region (${\mathrm{\Omega }}_1$) and the second penalizes distance between the desired minimizer $\mathbf{\Psi }$ and some negative μ inside ${\mathrm{\Omega }}_2$. A desired smoothness of the resulting level set function is obtained by searching for the minimizer inside a suitably chosen Sobolev space $W_{\alpha ,\beta }$, as introduced in Section 3.

We minimize the cost functional in (19(19) ${\varphi }_{new}=\underset{\mathrm{\Psi }}{\arg \min }K(\mathrm{\Psi })=\frac{{\alpha }_1}{2}||(\mathrm{\Psi }-{\varphi }_{old}){\chi }_1|{|}^2+\frac{{\alpha }_2}{2}||(\mathrm{\Psi }-\mu ){\chi }_2|{|}^2,$(19) ) for randomly selected locations for ${\mathrm{\Omega }}_2$ at different steps of our colour level set inversion routine, but only during a specific range of iteration numbers which defines the seeding phase. This process is external to the actual colour level set inversion routine and can be seen as off-line iterations which are performed either before or after selected updates for $\mathbf{\Phi }$. In context of the level set functions themselves, all we want to achieve in those external iterations is to approximate a given value by the level set function inside the seed region ${\mathrm{\Omega }}_2$, with minimal changes outside of it.

We summarize the new shape-based reconstruction approach discussed so far in form of a pseudo-code in Algorithm 1.

6. Numerical experiments

In this section, we will demonstrate the performance of this novel algorithm in a number of different situations, and compare it with a selection of alternative algorithms. In particular, we want to use the LK-colour level set scheme in Algorithm 1 for the application of imaging shielded regions. We choose to model the entire experimental domain as a cube-shaped region ${\mathrm{\Omega }}^h=[3]\times [3]\times [3]{\mathrm{m}}^3$ which includes a shielded container of size $[3]\times [3]\times [2.4]m^3$, similar to the situation depicted in Figure 2. As indicated in that figure, the regions of heights 0.3m above and below the container are reserved to accommodate the (source and receiver) antennas. The walls of the container are 0.3m thick at each side. We discretize this experimental domain ${\mathrm{\Omega }}^h$ by using a regular $20\times 20\times 40$ rectangular grid, each grid cell having physical dimension $[0.15]\times [0.15]\times [0.075]{\mathrm{m}}^3$. This means that the actual container including shields occupies $20\times 20\times 32$ grid cells. The internal region of the container without the shields is represented by $16\times 16\times 24$ grid cells, which amounts to a total of 6144 grid cells. Those are the locations where the physical parameter (conductivity) needs to be specified by the algorithm. In order to do so, we consider a number of $n_s=16$ sources and $n_r=16$ receivers, with all sources and receivers representing wire loops as indicated as well in Figure 2. Each source and receiver has dimension $[0.6]\times [0.6]\times [0.075]m^3$, with the sources being excited with an electric current $I_j=I=0.1A$ and each receiver recording a data value represented by (3(3) ${\left[{\mathcal{M}}_j{\mathbf{E}}_j\right]}_k={\oint }_{S_k}{\mathbf{E}}_j\cdot \hat{\mathbf{t}}\mathrm{d}l,$(3) ). A fixed and given probing frequency of f = 1MHz is applied in all cases.

Figure 2. Opposing arrays of sources and receivers with shielded cage in between. The interior of the box is hidden behind the shields. The sources and receivers are only located on top and bottom of the box, but other setups are also possible. Left: Complete set up of imaging the contents of a shielded cage. Right: Schematic of shielding in z direction. The domain of interest is between these shields in each direction.

Figure 3. 2D Visualization of level set representation described in (10(10) $b(\mathrm{\Phi })(\mathbf{x})=\{\begin{array}{ll}{b}^1& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})>0,\end{array}$(10) ); In this schematic visualization, it is assumed that the individual regions $\{\mathbf{x}:{\varphi }_1(\mathbf{x})\le 0\}$ and $\{\mathbf{x}:{\varphi }_2(\mathbf{x})\le 0\}$ are both disc-shaped, but that the specific rule (10(10) $b(\mathrm{\Phi })(\mathbf{x})=\{\begin{array}{ll}{b}^1& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1(\mathbf{x})>0\mathrm{and}{\varphi }_2(\mathbf{x})>0,\end{array}$(10) ) specifies which value to take in the overlap region $\{\mathbf{x}:{\varphi }_1(\mathbf{x})\le 0\text{ and }{\varphi }_2(\mathbf{x})\le 0\}$

On a more technical note, the computational domain is furthermore surrounded in our computational setup by additional absorbing boundary layers in order to prevent any outgoing fields from being reflected back into the computational domain. Our discretization of Maxwell's equations in (1a(1a) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{E}}_j(\mathbf{x})-a(\mathbf{x}){\mathbf{H}}_j(\mathbf{x})& ={\mathbf{M}}_j(\mathbf{x});\end{array}$(1a) ), (1b(1b) $\begin{array}{rl}\mathrm{\nabla }\times {\mathbf{H}}_j(\mathbf{x})-b(\mathbf{x}){\mathbf{E}}_j(\mathbf{x})& ={\mathbf{J}}_j(\mathbf{x}),\end{array}$(1b) ) follows closely the scheme described in [63–65] and has been implemented for our purposes in Python. It is a finite volume method using a Yee grid which is well suited to deal with high-contrast situations as the ones considered here. It solves for the EM fields $\mathbf{E}$ and $\mathbf{H}$ indirectly via finding a vector potential solution. The discrete Maxwell system is iteratively solved using the BiCGSTAB algorithm, following advice from [66], with a tolerance ${10}^{-4}$. For more details on the numerical code for solving Maxwell's equations in our numerical experiments, we refer to [48].

Here we demonstrate the performance of Algorithm 1 by using two numerical experiments, both of which involve placing two objects with differing constant conductivities in a background conductivity profile resembling air. The first is an object embedded inside another, whereas the second consists of two isolated objects. Both examples come with particular challenges.

In all numerical experiments, we choose this domain of interest to be shielded by a closed solid cage with conductivity of the walls being $b=0.1$ Sm${}^{-1}$. Currently, a level set representation is not used for the shields since they are assumed to be given and known. These are added to the conductivity profile once it has been defined by ${\varphi }_1$ and ${\varphi }_2$, for both data generation and when the level set functions are updated. As a default, we initialize both level set functions in the vector ${\mathrm{\Phi }}^{(0)}$ as ellipsoids in the central area of our domain (but not identical). This initial guess is by far not optimal and better initializations can be designed by simple preprocessing steps. However, for our purpose, an almost centred initial guess for each level set function suffices for demonstrating the general behaviour of the algorithm. Specific values for the conductivities used in the two numerical experiments are shown in Table 2.

Table 2. Conductivity parameters used in numerical experiments.

	$b^1$ (Sm${}^{-1}$)	$b^2$ (Sm${}^{-1}$)	$b^3$ (Sm${}^{-1}$)
Numerical Experiment 1	0.5	10.0	$1\times {10}^{-8}$
Numerical Experiment 2	0.5	1.0	$1\times {10}^{-8}$

In all our numerical experiments presented in this article, we choose to use synthetically generated data using the same forward modelling method, with additional noise added in the form of white Gaussian noise (1%). Therefore, the data ${\tilde{\mathcal{F}}}_j$ admits a decomposition (20) ${\tilde{\mathcal{F}}}_j={\mathcal{M}}_j{\hat{\mathbf{E}}}_j(\mathbf{x})(1+c{ϵ}_j),$(20) where ${\hat{\mathbf{E}}}_j(\mathbf{x})$ is the electric field generated by the true objects in each numerical experiment, ${\mathcal{M}}_j$ is given by (3(3) ${\left[{\mathcal{M}}_j{\mathbf{E}}_j\right]}_k={\oint }_{S_k}{\mathbf{E}}_j\cdot \hat{\mathbf{t}}\mathrm{d}l,$(3) ), ${ϵ}_j\sim \mathcal{N}(0,1)$ and $c\in \mathbb{R}$ is a scaling parameter representing the noise level. In both numerical experiments we assume the computational setup described above and that the measured data ${\tilde{\mathcal{F}}}_j$ used for recovery follows the form described in (20(20) ${\tilde{\mathcal{F}}}_j={\mathcal{M}}_j{\hat{\mathbf{E}}}_j(\mathbf{x})(1+c{ϵ}_j),$(20) ). Notice also that, in order to avoid the so-called ‘inverse crime’, we have included three different variations of the first numerical experiment where one of the background distributions occupying one of the regions is perturbed by different types of additional (including non-Gaussian) noise. We will comment on that further below. In our future work, we plan to use data generated by different grid sizes and/or different forward solvers, but for our proof-of-concept study presented here, we are confident that adding Gaussian noise fits the same purpose and allows us to evaluate better the quality of the final reconstructions against more traditional inversion schemes.

6.1. LK-colour level set inversion for numerical experiment 1

We refer to Appendix for details on the particular values chosen for the line search scheme. Even though the line search criteria are quite technical, they form an important part of our reconstruction algorithm. Therefore, along with images of the colour level set reconstruction, we also display information on the step size in relationship with the average number of voxels that change per sweep and information on how the pseudo-cost behaves (the concept of a pseudo-cost is defined further below when discussing Figure 5 in more detail). In the corresponding Figure 5 displaying the step size information, we display the admissible interval range ${\hat{I}}_{target}^m$ as solid straight lines (with black colour in the online version of this paper), and bounds of the desired region ${\hat{I}}_{target}^{\prime m}$ are shown in dashed straight lines (with green colour in the online version).

Figure 4 shows various stages of a 3D LK-colour level set inversion for the Numerical Experiment 1 in form of surface plots. The bottom right image displays the true phantom, and the upper left image shows the initial guess. The remaining images display snapshots of the shape evolution at various sweep numbers (each sweep considering all source positions exactly once). Cross-sections of the final result are also displayed in the fourth row of Figure 6 which will be discussed further below. It becomes visible that, upon convergence, the algorithm has managed to capture the main topological characteristics of the true phantom, even though we have by far less data available than we have pixels with unknown conductivity values. This is even more impressive considering that each of the three conductivity regions present in the true phantom has quite a high contrast with respect to the others.

Figure 4. Surface plots of 3D shape evolution for Numerical Experiment 1. Shown are the initial shape and snapshots at iteration numbers $s=5,25,50,100,200,400,450$ and the true phantom. The surrounding shields are not shown here. In the online version of this paper, the  colour red indicates the shape of the conductivity $b^2$ and blue indicates the shape of the conductivity $b^1$. The third region is transparent and represents air ($b^3$).

Figure 5. Line search analysis of LK-colour level set reconstruction for Numerical Experiment 1. Top row: Evolution of the average number of voxels $\overline{N}$ as a function of the sweep number s. The dotted (and yellow in the online version of the paper)  line represents ${\overline{N}}_s^1$ and the solid (black in the online version)  line represents ${\overline{N}}_s^2$. Middle row: Evolution of the step sizes ${\tau }_1^s$ and ${\tau }_2^s$ , denoted in the online version by yellow and black lines. Bottom row: Evolution of pseudo-cost $\tilde{J}[\mathrm{\Phi }{]}^s$ as a function of the sweep number s.

Figure 6. 2D cross-sections through a 3D reconstruction for Numerical Experiment 1 using different regularization schemes. Shown is the container region including shields. Each row: Left z = 25, middle $x=15$, right y = 18. 1st row: LK-Pixel, 2nd row: LK-single-level set $b_i=\frac{1}{2}(b^1+b^2)$, 3rd row: LK-single-level set $b_i=b^2$, 4th row: LK-colour level set and 5th row: true phantom.

Compare already here the corresponding evolution of line search criteria and cost functional displayed in Figure 5 (discussed further below). As expected, the reconstruction cannot be perfect even though the value of the cost functional has been reduced significantly by sweep s = 500. Based on the above count of unknowns and data, one expects some form of non-uniqueness in the reconstructions which cannot be removed without additional prior knowledge. Indeed, the final reconstruction shows a higher amount of fine details compared to the true phantom. Some smaller parts of value $b^2$ are attached to the main body of value $b^1$. This appears to be an artefact that stems from details of the shape evolution and can be characterized as a typical local minimum. Removing those artefacts would not make a significant change in the final cost. On the other hand, it is comforting to observe that the annulus in the region of $b^2$ is correctly identified, even though it is not included in the initial guess. Overall most of the material with values $b^1$ and $b^2$ are contained in the correct area of the container when compared with the true model.

Figure 5 shows reconstruction information of Numerical Experiment 1. In the top figure, we observe the average number of voxels that change each sweep for both level set functions (dotted line (yellow in the online version): ${\varphi }_1$, solid line (black in the online version): ${\varphi }_2$). The top solid straight line descends as the algorithm progresses, making the admissible interval smaller. This induces ${\hat{I}}_{target}^{\prime m}$ to also be smaller, meaning acceptance of an update is more difficult.

Rather than monitoring the actual cost, which is computationally demanding, we choose to track a different measure which is cheap to compute. As indicated in Algorithm 1, let s denote the sweep number, j the index of the source position and $J_j[{\mathrm{\Phi }}_j^s]$ the corresponding cost value as calculated according to (12(12) $\mathcal{T}(\mathrm{\Phi })=\mathcal{R}[b(\mathrm{\Phi })],\mathcal{J}(\mathrm{\Phi })=||\mathcal{T}(\mathrm{\Phi })|{|}_{{\mathcal{L}}^2(\mathrm{\Omega })}^2,$(12) ) in order to obtain the next update. The pseudo-cost ${\tilde{J}}_s$ and actual cost ${\hat{J}}_s$, both obtained upon completion of sweep number s, are defined as (21) ${\tilde{J}}_s=\sum _{j=0}^{n_s-1}{J}_j[{\mathrm{\Phi }}_j^s],{\hat{J}}_s=\sum _{j=0}^{n_s-1}{J}_j[{\mathrm{\Phi }}_{(n_s-1)}^s]=J\left({\mathrm{\Phi }}_{(n_s-1)}^s\right).$(21) The actual cost ${\hat{J}}_s$ corresponds to the standard cost value evaluated off-line by using the last level set function obtained by sweep s. Its calculation requires us to run $n_s$ forward problems using our Maxwell solver on this single-level set function. The pseudo-cost ${\tilde{J}}_s$, on the other hand, uses partial cost values associated with the level set functions that are calculated by Algorithm 1 during sweep s in order to obtain the individual updates. Its calculation comes without any significant additional expenses. It can be viewed as a lagged cost and is significantly less expensive to obtain than computing the actual cost ${\hat{J}}_s$ after each sweep. When adopting our line search strategy, significant differences between the pseudo and actual cost have not been observed throughout the evolution when choosing to calculate both in the numerical experiments. Therefore, only ${\tilde{J}}_s$ is calculated and displayed here.

In this particular example, now observing the middle figure, we see that the algorithm has managed to find a step size which allows ${\overline{N}}^m\in {\hat{I}}_{target}^{\prime m}$ for the majority of sweeps. The bottom figure shows the pseudo-cost in each sweep, which decreases as the algorithm progresses, meaning the error between the data generated by our reconstruction and the test phantom is becoming smaller.

6.1.1. Comparison of regularization schemes

Figure 6 shows a comparison of the colour level set inversion scheme with other reconstruction algorithms for Numerical Experiment 1. As part of this, we compare our colour level set reconstructions with two variants of a single-level set inversion and also a traditional ${\mathcal{L}}^2$-pixel-based scheme. The two single-level set inversions differ in their a priori information; the first (diplayed in the second row of the figure) assumes that the single interior conductivity value in the single-level set inversion, which we label $b_i$, is an arithmetic average of the two different values present in the true phantom. The second (displayed in row three of the figure) assumes that it coincides with the highest of the two values, namely $b_i=b^2$. The exterior background in both these cases is assumed to be $b^3$, as in the colour level set approach.

The figure moreover shows in row one that the ${\mathcal{L}}^2$-pixel-based scheme performs reasonably well given that the true phantom has relatively high contrast between conductivity domains. This reconstruction roughly corresponds to the upper left panel of the flow chart provided in Figure 1. It satisfies the data but lacks the clear structure suggested by our a priori information. It can be argued whether applying a standard off-the-shelf image segmentation approach to the standard Tikhonov Philips type result shown in row one of this Figure 6 would be a good approximation to the correct reference figure, which is shown in row five of the same figure. As mentioned before, the pixel-based scheme promotes over-smooth conductivity profiles in low-frequency situations, whereas in contrast the level set inversion schemes enforce non-smooth conductivity profiles with sharp edges. When given correct a priori information on conductivity values, the level set inversion can resolve interfaces between shapes more clearly.

However, in this more challenging case, the binary approach using just one level set function clearly might not be the best possible approach either, as seen in rows two and three of the figure. Apparently, more complicated scenarios such as the true phantom for Numerical Experiment 1 do not lend themselves well to the single-level set inversion, since fundamentally the method is designed to only recover shapes of two-parameter domains. This observation motivates the extended studies performed in this paper on employing more general colour level set schemes instead.

But before examining results of the novel colour level set approach, let us compare two values for $b_i$ in this a priori incorrectly designed level set model. As shown in [4], using an incorrect value for the level set inversion tends to deliver growing or shrinking domains to fit the data depending on whether a priori information on the conductivity is lower or higher than what is present in the true phantom. We observe this phenomenon in the LK-single-level set reconstruction here. When we take the interior conductivity value to be an average of the two present in the true phantom, we observe an inflated or deflated conductivity profile as the algorithm tries to adjust the shape to fit the data. We can compare this situation to the second scenario, where we take the interior conductivity value to be the highest interior value present in the true phantom. With this choice, we obtain a much smaller interior shape since the interior conductivity is almost double that of the averaged interior value.

The results of the LK-colour level set inversion algorithm, shown in row four of the Figure 6, demonstrate that this novel scheme performs much better than the standard level set approach (and also better than the pixel-based approach) in attempting to recover the true phantom. This is so since characteristics of the true phantom are present in the reconstruction. For example, the scheme has managed to recover a conductivity profile that has the same distinctive separation features (i.e. a portion of $b^3$ is trapped in between $b^1$ and $b^2$). In comparison, the pixel-based scheme delivers a smooth conductivity profile extended over the entire domain, and the single-level set scheme delivers an inflated or deflated shape (depending on the single internal value used) of the true phantom.

Figure 7 shows a comparison of the pseudo-cost for each reconstruction scheme shown in Figure 6. The ${\mathcal{L}}^2$-pixel-based scheme has the highest data misfit value, whereas the data misfit of both single-level set schemes are relatively close to the colour level set scheme despite poor recovery of the true phantom. The relatively high value for the pixel-based scheme is a little bit surprising, but can be due to very slow convergence of this scheme. Displayed is the value after s = 500 sweeps, but it is possible that further reduction is achieved by keeping the scheme running significantly longer (which however would be impractical in realistic applications). It might also be due to the smoothing effect of the pixel-based reconstruction schemes, which makes it difficult to correctly fit data obtained with a high-contrast model from reconstructions of very low contrast. The two-value level set approach performs better here, but still cannot match the low final cost value of the colour level set reconstructions in this example. This is certainly plausible, since the (three-value) colour level set model agrees best of all with the correct setup of the true phantom and incorporates optimally the available prior information.

Figure 7. Pseudo-cost for four different regularization schemes; LK-colour level set, LK-single-level set, LK-pixel and SD-pixel.

6.1.2. Variability in the true conductivity profile

We can test the numerical stability of the colour level set scheme by increasing the realism of the true phantom inside the container. In this case, we assume that the background conductivity $b^1$ (other than air which is $b^3$) has some variability. Here, we assume that the true conductivity profile admits the following decomposition: (22) $b({\mathrm{\Phi }}^{\ast })(\mathbf{x})=\{\begin{array}{ll}\xi (\mathbf{x})& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})>0,\end{array}$(22) where $\xi (\mathbf{x})$ is drawn from a probability distribution and the ${}^{\prime }{\ast }^{\prime }$ superscript denotes the (discretized) true level set function. In other words, the parameter content of region $S_1$ is actually random with properly chosen statistical characteristics, to be explained below. Those distributions are chosen to depict possible variations of actual boxes or cargo containers in realistic situations. Here we will consider a comparison between two variations of a normal distribution N and a uniform distribution U for the true content of region $S_1$. Therefore, in the two cases, $\xi ({\mathbf{x}}_i)$ admits the decompositions (23a) $\begin{array}{rl}\xi ({\mathbf{x}}_i)& \sim N(\beta b^1,\alpha b^1)\end{array}$(23a) (23b) $\begin{array}{rl}\xi ({\mathbf{x}}_i)& \sim U(b^1-\delta ,b^1+\delta ),\end{array}$(23b) where $\alpha ,\beta ,\delta \in {\mathbb{R}}^{+}$ are chosen parameters. Note that our two variants of drawing from a normal distribution in (23(23a) $\begin{array}{rl}\xi ({\mathbf{x}}_i)& \sim N(\beta b^1,\alpha b^1)\end{array}$(23a) ) stem from choosing various α which determines the variance.

Numerically speaking we now assign each cell in $S^1$ of the true phantom to a single draw from a probability distribution without replacement. We consider three examples; the first where $\xi (\mathbf{x})$ is drawn from a normal distribution such that the conductivity associated with $S^1$ is now slowly varying around $b^1$; the second where $\xi (\mathbf{x}\mathbf{)}$ is drawn from a uniform distribution whose mean is $b^1$; and the third where $\xi (\mathbf{x})$ is drawn from a normal distribution with mean $2b^1$. See Figure 8.

Figure 8. Left: Sample of conductivity profile in $S^1$ using probability distribution $N(b^1,0.05)$, Middle: Sample of conductivity profile in $S^1$ using probability distribution $U(0.2,0.8)$, Right: Sample of conductivity profile in $S^1$ using probability distribution $N(2b^1,0.25)$.

Figure 8 shows the statistical samples generated for building the true conductivity profiles associated with $S^1$ used in the inversions depicted in Figures 9– 11, from left to right, respectively. In all three figures, the  horizontal (red in the online version) line denotes the value of conductivity $b^1$ which is used in the colour level set inversion scheme.

Figure 9. 2D cross-sections through a 3D reconstruction for Numerical Experiment 1 using the formulation in (22(22) $b({\mathrm{\Phi }}^{\ast })(\mathbf{x})=\{\begin{array}{ll}\xi (\mathbf{x})& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})>0,\end{array}$(22) ) for the true phantom. Here, $\xi (\mathbf{x})$ is drawn from a normal distribution with parameters $\beta =1,\alpha =10$. Each row: Left z = 23, middle x = 15, right y = 18 Top: true phantom, Bottom: LK-colour level set reconstruction at sweep s = 1000.

Figure 10. 2D cross-sections through a 3D reconstruction for Numerical Experiment 1 using the formulation in (22(22) $b({\mathrm{\Phi }}^{\ast })(\mathbf{x})=\{\begin{array}{ll}\xi (\mathbf{x})& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})>0,\end{array}$(22) ) for the true phantom. Here, $\xi (\mathbf{x})$ is drawn from a uniform distribution with parameter $\delta =0.3$. Each row: Left z = 23, middle x = 15, right $y=18$ Top: True Phantom, Bottom: LK-colour level set reconstruction at sweep s = 10000

Figure 11. 2D cross-sections through a 3D reconstruction for Numerical Experiment 1 using the formulation in (22(22) $b({\mathrm{\Phi }}^{\ast })(\mathbf{x})=\{\begin{array}{ll}\xi (\mathbf{x})& \mathrm{in}{S}_1\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})\le 0,\\ b^2& \mathrm{in}{S}_2\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})\le 0,\\ b^3& \mathrm{in}{S}_3=\mathrm{\Omega }\mathrm{\setminus }(S_1\cup S_2)\mathrm{where}{\varphi }_1^{\ast }(\mathbf{x})>0\mathrm{and}{\varphi }_2^{\ast }(\mathbf{x})>0,\end{array}$(22) ) for the true phantom. Here, $\xi (\mathbf{x})$ is drawn from a normal distribution with parameters $\beta =2,\alpha =2$. Each row: Left z = 23, middle x = 15, right y = 18 Top: true phantom, Bottom: LK-colour level set reconstruction at sweep s = 1000.

Let us now apply the colour level set inversion scheme to this new challenging situation where content in the region $S_1$ is still assumed constant during the reconstruction, reflecting our lack of knowledge of the detailed nature of the content inside the box. Figures 9–11 depict LK-colour level set reconstructions where the true $S^1$ is characterized by using draws from both normal and uniform distributions.

All three figures show a LK-colour level set reconstruction where the region $S^1$ in the true phantom has a significant variability of internal conductivity associated with it. In comparison to the reconstructions shown in both Figures 4 and 6, we see that the reconstruction of $S^3$ is no longer occupying significant parts of the central region between the structures described by $S^1$ and $S^2$ (representing a gap filled with air in the true phantom). However, $S^2$ still has resemblance with the true phantom and some characteristics remain from the reconstruction involving a constant conductivity profile associated with $S^2$.

Apparently the samples shown in the left and middle images of Figure 8 do not inhibit the LK-colour level set inversion in its recovery of the shape in Figures 9 and 10 since values in $S^1$ fall either side of $b^1$. One possible explanation of this behaviour is that the underestimation and overestimation of the conductivity associated with $S^1$ in the colour level set inversion effectively cancel each other out when observing the data (the growing and shrinking effect as discussed earlier), since both these reconstructions are very similar to that shown in Figure 4.

In contrast, the conductivity reconstruction in Figure 11 associated with $S^1$ for the level set inversion using the sample in the right image of Figure 8 does not perform as well as the others. Even though the reconstruction overall captures the true behaviour of the conductivity profile, its resemblance to the true phantom is slightly worse than the other two sampling techniques, as it should be expected. This most likely occurs because the a priori information given to the algorithm overall underestimates conductivity values in large parts of the domain. Though it is encouraging that even in penalized situations like this the algorithm can still capture many characteristics which are truly present.

Nevertheless, all three reconstructions are impressive considering many of the conductivity values associated with $S^1$ are not used as a priori information in the LK-colour level set reconstruction. Figure 11, for example, shows a LK-colour level set reconstruction of a true phantom whereby the conductivity profile associated with $S^1$ has large variability and has mean around twice the true value. Although the reconstruction scheme has been heavily penalized with wrong a priori information, the recovery is quite similar to the first two, which is impressive.

Figure 12 shows a comparison of the pseudo-cost for each sampling technique.

Figure 12. Pseudo-cost comparison between different sampling methods for creating conductivity profile associated with $S^1$.

6.2. LK-colour level set inversion for numerical experiment 2

After the very convincing results seen in Numerical Experiment 1, we want to present now a different setup which highlights some potential pitfalls of the newly proposed algorithm. We also present a way of circumventing those maybe unexpected difficulties. In a similar way as for Numerical Experiment 1, Figure 13 shows various stages of a 3D LK-colour level set reconstruction for Numerical Experiment 2. The general setup of this example is similar to the one in Numerical Experiment 1, but the true profile (shown in the bottom right image of Figure 13) is quite different here with both objects being separated by a very large distance.

Figure 13. Surface plots of 3D shape evolution for Numerical Experiment 2. Shown are the initial shape and snapshots at iteration numbers $s=5,25,200,500$ and the true phantom. The surrounding shields are not shown here. In the online version of this paper, the colour red indicates the shape of the conductivity $b^2$ and blue indicates the shape of the conductivity $b^1$.

By observation, the algorithm performs well in the sense of locating both objects inside the imaging domain. However, the classification of the objects is obviously biased by choice of the initial guess and is done incorrectly in this particular reconstruction. This is problematic in situations where the exact nature of the recovered objects is of importance. Notice that this incorrect classification is compensated by the algorithm because it provides inflated or deflated areas for each of the two conductivity domains as a result of optimally fitting the cost functional.

We mention here without showing this result (for the sake of brevity), that if the initial guess is flipped, in the sense that ${\varphi }_1^{(0)}$ was changed to ${\varphi }_2^{(0)}$ and vice-versa, we recover the correct classification and location as it might be expected. Accordingly, choice of the initial guess can significantly influence the final reconstruction in the colour level set approach, even more than it is usually seen in the classical two-value level set method. This indicates that expectations on distinguishing between all regions correctly should not be set too high in this approach. The reconstructed internal values at each location can be affected by a new type of local minima in the chosen optimization approach, not present in the other methods discussed in this paper. We emphasize, though, that reconstructions using the colour level set scheme still seem at least of the same quality and usefulness as the more traditional methods compared to in this paper. None of those would be able to correctly identify the actual contrast to the background of each object.

Certainly, some of the previously emphasized additional advantages of our novel scheme are not guaranteed here and depend on choice of the initial guess. Certainly, remedies are available and we will present one suggestion using random seeding of objects next. Alternatively, more sophisticated global optimization approaches (as for example ensemble-based methods [61]) might also yield more reliable results for the colour level set scheme, for the price of increased computational expenses. Those are left for future research.

The proposed seeding method integrated into the colour level set scheme represents a kind of ‘globalization’ of the local descent technique described so far. It adds a seeding stage into the algorithm which can be interpreted as providing small perturbations or ‘jumps’ across the profile of the cost functional which can help guiding the scheme out of a local minimum. In more detail, we choose throughout a certain stage of the inversion process small parts of the imaging domain to be randomly selected for placement of an artificial seed object at different iterations. In our numerical experiment for testing this seeding technique, we choose each artificial object to have length 6 cells in each direction, resembling a cube seed. This represents $3\mathrm{\%}$ of the imaging domain. For comparison with Figure 13, we use the same initial guess. Figure 14 shows various stages of a 3D LK-colour level set reconstruction with an added seeding phase, where $\gamma =15$, but keeping the misleading initial starting model. The reconstruction is an improvement on that in Figure 13, in both classification and location. Therefore, the added seeding phase has been rewarded by reaching a ‘better local minimum’ (closer to the desired ‘global minimum’) judged by visual inspection when knowing the true phantom.

Figure 14. Surface plots of 3D shape evolution for Numerical Experiment 2 with an initial seeding phase where $\gamma =15$. The surrounding shields are not shown here. Displayed are the initial shape and snapshots at iteration numbers s = 5, 25, 200, 500 and the true phantom. In the online version of this paper, the colour red indicates the shape of the conductivity $b^2$ and blue indicates the shape of the conductivity $b^1$.

Certainly, in practical scenarios, the true phantom is not known and a more systematic analysis of the benefits of such a seeding phase would be desirable. In order to understand those benefits better, we have performed many simulations to see and compare their rewards in a statistical setting. Running 125 simulations of the inversion with seeding (still keeping the biased initial starting model) revealed that only $15\mathrm{\%}$ of reconstructions gave correct classification of objects. In another run of numerical experiments without a biased initial guess, where the initialization would completely rely on seeding objects throughout the reconstruction, we have seen an approximately 50/50 chance of reconstructing the correct shape, as it would be expected.

For this particular experiment of two isolated objects of small to medium size, we have not found a clear measure for indicating in the final result which one represents the correct contrast values for each of them. The main difficulty, in addition to obviously representing local minima of the underlying cost functional, is that the final cost value for both reconstructions does not show significant differences between those different reconstruction options. This is due to the mutual compensation of the object's assumed internal value and its reconstructed size via the already mentioned inflation and deflation mechanism in the inversion process. This behaviour reflects the severe ill-posedness of the underlying inverse problem. It does affect certain geometrical setups of the true phantom more than others, as seen in the first numerical experiment where this did not cause any significant difficulties.

Even though the presented seeding phase does not always guarantee to reach the global minimum of the chosen cost functional, as demonstrated in this particular numerical example, we want to emphasize that it usually does improve results and is useful for speeding up convergence in most cases. Moreover, it is quite flexible and can be administered at any stage of the inversion routine. In [6], they decide to seed in various cycles. We decide to have an initial burn-in phase, whereby we seed both ${\varphi }_1$ and ${\varphi }_2$ at the beginning of each sweep s for $0\le s<\gamma <S$. This amounts to solving γ small-scale minimization problems of (19(19) ${\varphi }_{new}=\underset{\mathrm{\Psi }}{\arg \min }K(\mathrm{\Psi })=\frac{{\alpha }_1}{2}||(\mathrm{\Psi }-{\varphi }_{old}){\chi }_1|{|}^2+\frac{{\alpha }_2}{2}||(\mathrm{\Psi }-\mu ){\chi }_2|{|}^2,$(19) ) in the burn-in phase. Choosing this specific setup arises from an observation: when level set functions have converged in the colour level set regime, it tends to be difficult for the shapes to be flipped in value by seeding. Therefore, allowing the process to be early in the inversion is more favourable if we are to overcome this difficulty. Moreover, the seeding phase can be viewed as a dynamic initial guess scheme to encourage the minimization process to leave the vicinity of a local minimizer resembling incorrect classification. In this sense, it can be considered a ‘low cost version of globalization’ of the applied local optimization scheme.

Employing more sophisticated truly global optimization schemes on colour level set models is possible as well, but we expect them to be computationally far more expensive. As part of such approaches, ensemble-based approaches would be a natural option where several possible final results are recorded and statistically evaluated [61]. We also would like to mention that in many practical application it is often not really necessary to identify one particular global minimum but rather to be able to classify a small number of possible local minima against certain criteria, e.g. identifying a possible shape hidden in a box or in a cargo container. Since the shapes of the expected local minima in our examples are all very similar, such a classification can be performed reliably in most cases.

Finally, also a joint reconstruction of all three internal parameter values and the corresponding shapes could be attempted in order to address this challenging situation. This is in principle possible by following the general lines of theory as demonstrated for example in [18,23,67]. However, this process of joint shape and parameter inversion tends to be extremely slow in this low frequency situation, and appropriate line search strategies for efficiently combining the searches of those quite different quantities (shape and contrast value) in our situation are still missing.

7. Summary and conclusions

We have introduced a novel regularization scheme for imaging shielded regions filled with high contrast conductivity distributions using low frequency electromagnetics in 3D. The LK-colour level set scheme has been shown to produce interesting and in many situations significantly improved results compared to more classical approaches. We provided a comparison with alternative regularization schemes for this scenario to see how well it fared against existing methods. The scenario considered in Numerical Experiment 1 is more realistic and challenging than those previously considered for this application. Increasing realism of the problem by choosing a conductivity profile with more than two constant domains as the true phantom, has led to improved results but also to additional challenges. In the LK-colour level set regime, we found that the increased complexity of the model comes with some new types of local minima. Some of those local minima resemble incorrect classification of internal conductivity values of each individual object compared to the true value, regardless of a good approximation of its shape. Here, the seeding process can try to circumvent some of those issues, or at least give a set of local minima to infer from by other means. Alternative ensemble-based methods might provide a better global approach for characterizing such situations at the expense of higher computational cost.

For future research, it has been shown that in most cases the use of more than one level set function clearly demonstrates great promise when applied to a variety of particular situations. We have used this extended flexibility here for modelling more than two regions of different internal conductivity values. In addition to this, it can also be applied to a joint inversion for more than one physical parameter. Clearly, different level set functions can not only represent different values of one given physical parameter, but in addition they can represent different physical parameters. For example, another material parameter variation such as the magnetic permeability μ occupying a different region other than the conductivity, can be incorporated into the model for the reconstruction. In particular, for the application of screening boxes or containers, this might be of practical significance. Anisotropic behaviour of some EM parameters might occur in realistic situations, which in principle can as well be addressed by a similar algorithm as the one presented here. Also, properties of the shields can be incorporated into the inversion scheme and the use of data from more than one frequency can be employed. Finally, the combined reconstruction-segmentation approach presented here is not limited to low-frequency electromagnetic imaging. It can be applied to a large variety of different imaging technologies where it might provide very promising advantages compared to the more traditional approach of first reconstructing and then segmenting the obtained results.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The work of A. J. H. was funded by the UK Ministry of Defence, as part of his Ph.D. project. Ministry of Defence © British Crown Owned Copyright 2020/AWE.

Appendix. LK-colour level set line search criteria

In this Appendix, we introduce a line search technique for a colour level set inversion involving N level set functions. The scheme is used to choose a step size vector ${\mathit{\tau }}^s=[{\tau }_1^s,{\tau }_2^s,\dots ,{\tau }_N^s]$ in the update formula (A1) ${\mathrm{\Phi }}_{j+1}^s={\mathrm{\Phi }}_j^s-{\mathit{\tau }}^s\circ \mathrm{\nabla }{\mathcal{J}}_j[{\mathrm{\Phi }}_j^s],$(A1) where s denotes the sweep number.

Choosing the step size vector ${\mathit{\tau }}^s$ is tricky. Here, we adopt ideas from [4] and choose to modify (A1(A1) ${\mathrm{\Phi }}_{j+1}^s={\mathrm{\Phi }}_j^s-{\mathit{\tau }}^s\circ \mathrm{\nabla }{\mathcal{J}}_j[{\mathrm{\Phi }}_j^s],$(A1) ) such that (A2) ${\mathit{\tau }}^s=\mathit{f}({\mathit{w}}^s\circ {\mathit{\tau }}^0),$(A2) where $\mathit{f}(\cdot )$ is the vector-valued form of the backtracking scheme given in Algorithm 2, ${\mathit{w}}^s=[w_1,w_2,\dots ,w_N]$ is a vector of scalars which are adjusted at the end of each sweep and ${\mathit{\tau }}^0$ is a vector of initializations of the step size for each level set function in $\mathbf{\Phi }$. Hence, we arrive at the update formula: (A3) ${\mathrm{\Phi }}_{j+1}^s={\mathrm{\Phi }}_j^s-\mathit{f}({\mathit{w}}^s\circ {\mathit{\tau }}^0)\circ \mathrm{\nabla }{\mathcal{J}}_j[{\mathrm{\Phi }}_j^s].$(A3)

The step size in (A3(A3) ${\mathrm{\Phi }}_{j+1}^s={\mathrm{\Phi }}_j^s-\mathit{f}({\mathit{w}}^s\circ {\mathit{\tau }}^0)\circ \mathrm{\nabla }{\mathcal{J}}_j[{\mathrm{\Phi }}_j^s].$(A3) ) has two components: the first is applied every update ($\mathit{f}$) and the second is a dynamic scheme which alters the input of $\mathit{f}$ at the end of each sweep.

To begin, we discuss the first component of the scheme. Here, our goal is to choose each entry in ${\mathit{\tau }}^s$ such that $N_j^m$ (the number of voxels that change per jth source for the mth level set function, $j=0,1,\dots ,n_s-1$) is contained in an interval which is deemed suitable for the shape associated with the mth level set function to evolve smoothly. The relationship between $N_j^m$ and the data misfit is of interest, but as of present is unclear. Although the relationship is unknown, intuition tells us that if we only allow a small amount of voxels to change per update, then we ensure a smooth evolution of the shape which also hopefully induces a smooth evolution of the cost. We choose $N_j^m\in I_{target}^m$, where $I_{target}^m=[N_{inf}^m,N_{sup}^m]$ is selected appropriately. Typically, $N_{inf}^m=0$ and $N_{sup}^m$ is chosen as a small percentage of the total number of voxels in the domain. This quantity gives us a metric for how much each level set function evolves per sweep.

Although the backtracking line search makes sure that updates are only made if $N_j^m\in I_{target}^m$, it does not ensure an unequal weighting to each sensitivity pattern generated by the jth source. Ideally, we would like to give greater weighting to those sources which capture objects present in their sensitivity patterns. For example, if we initialize the pth entry in ${\mathit{\tau }}^s$ such that an update to the pth level set function is large, we will consistently accept a voxel change just below $N_{sup}^p$ for all sources. This introduces two problems; the first is that more iterations of the backtracking scheme will have to be computed in order to reduce ${\tau }_p^0$ such that $N_{sup}^p\in I_{target}^p$ and the second is that we will give equal weighting to so-called bad updates in the Kaczmarz scheme.

To alleviate this problem, we make use of the quantity $N_j^m$. Once a sweep is complete (i.e. we make a Kaczmarz update for each source), we compute the average number of voxels that have changed for that sweep. Mathematically, we have that (A4) ${\overline{N}}^m=\frac{1}{n_s}\sum _{j=1}^{n_s}{N}_j^m.$(A4) This brings us to the second component of the line search criterion. We can adjust ${\mathit{\tau }}^0$ with a weighting vector ${\mathit{w}}^s$, which is chosen by some criterion involving (A4(A4) ${\overline{N}}^m=\frac{1}{n_s}\sum _{j=1}^{n_s}{N}_j^m.$(A4) ). Let us now introduce the subinterval of $I_{target}^{\prime m}=[N_{low}^m,N_{high}^m]$, which will be included in our criterion for updating the entries in ${\mathit{w}}^s$: (A5) $\begin{array}{rl}{N}_{high}^m& =(1-d_m)N_{sup}^m,\end{array}$(A5) (A6) $\begin{array}{rl}{N}_{low}^m& =d_m{N}_{sup},\end{array}$(A6) where $d_m\in (0,1/2)$. From this, we can deduce that (A7) $N_{inf}^m<N_{low}^m<N_{high}^m<N_{sup}^m.$(A7) Following this, we adjust each entry in ${\mathit{w}}^s$ according to some criterion involving (A4(A4) ${\overline{N}}^m=\frac{1}{n_s}\sum _{j=1}^{n_s}{N}_j^m.$(A4) ), (A5(A5) $\begin{array}{rl}{N}_{high}^m& =(1-d_m)N_{sup}^m,\end{array}$(A5) ), (A6(A6) $\begin{array}{rl}{N}_{low}^m& =d_m{N}_{sup},\end{array}$(A6) ). For brevity, we show how the pth component of ${\mathit{w}}^s$ is updated per sweep: (A8) $w_p^s=\{\begin{array}{ll}{k}_p^1{w}_p^{s-1}& \mathrm{if}{\overline{N}}_s^p>N_{high}^p\mathrm{where}{k}_p^1\in [0.5,1);\\ k_p^2{w}_p^{s-1}& \mathrm{if}{\overline{N}}_s^p<N_{low}^p\mathrm{where}{k}_p^2\in (1,2];\\ w_p^{s-1}& \mathrm{if}{N}_{low}^p\le {\overline{N}}_s^p\le N_{high}^p,\end{array}$(A8) where s denotes the quantity for sweep number s.

In other words, this means that if the pth entry to the backtracking line search in Algorithm 2 is too large then this criterion will reduce the step size before it is altered by the backtracking line search. Likewise, if the step size is too low, then this scheme can help the updates become larger. Of course, if entries are extremely large or small, then a scalar multiplier in the intervals given may be insufficient. Nonetheless, the scheme performs reasonably well for our application.

Note that so far, both intervals $I_{target}^m$ and $I_{target}^{\prime m}$ have constant lower and upper bounds. These can be generalized to include bounds which are linearly decreasing dependent on the sweep number. This is not arbitrary, as we expect diminishing returns of gradient descent methods as the sweep number increases. Insisting on the same amount of voxels to be altered later in the algorithm is likely to lead to artefacts in the reconstruction, as the data misfit reduces slower in later iterations. The proportion between $N_j^m$ and data misfit should remain the same, which is what these dynamic intervals hope to capture.

Therefore, we introduce a dynamic interval $I_{target}^m=I_{target}^m(s,{\eta }_m)$, which induces ${\hat{I}}_{target}^m={\hat{I}}_{target}^m(s,{\eta }_m)$. We choose the bounds of ${\hat{I}}_{target}^m$ to behave linearly as a function of the sweep number s. It follows that (A9) $\begin{array}{rl}{\hat{N}}_{low}^m(s,{\eta }_m)& =-\frac{1}{S}{N}_{low}^m{\eta }_{m}s+N_{low}^m,\end{array}$(A9) (A10) $\begin{array}{rl}{\hat{N}}_{high}^m(s,{\eta }_m)& =\frac{1}{S}(N_{low}^m-N_{high}^m){\eta }_{m}s+N_{high}^m;\end{array}$(A10) where ${\eta }_m\in [0,1]$. Defining the acceptance intervals ${\hat{I}}_{target}^m(s,{\eta }_m)$ and using them for the line search criterion is summarized below:

1. Choose $N_{inf}^m$ and $N_{sup}^m$.

2. Determine $N_{low}=c_m{N}_{sup}^m$ and $N_{high}^m=(1-c_m)N_{sup}^m$.

3. Compute ${\hat{N}}_{low}^m(s,{\eta }_m)$ and ${\hat{N}}_{high}^m(s,{\eta }_m)$ using (A9(A9) $\begin{array}{rl}{\hat{N}}_{low}^m(s,{\eta }_m)& =-\frac{1}{S}{N}_{low}^m{\eta }_{m}s+N_{low}^m,\end{array}$(A9) ), (A10(A10) $\begin{array}{rl}{\hat{N}}_{high}^m(s,{\eta }_m)& =\frac{1}{S}(N_{low}^m-N_{high}^m){\eta }_{m}s+N_{high}^m;\end{array}$(A10) ).

4. Determine ${\hat{I}}_{target}^m(s,{\eta }_m)$ using relationship in (2).

5. Use ${\hat{I}}_{target}^{\prime m}(s,{\eta }_m)=[{\hat{N}}_{low}^m,{\hat{N}}_{high}^m]$ to determine ${\mathit{w}}^s$ with criterion (A8(A8) $w_p^s=\{\begin{array}{ll}{k}_p^1{w}_p^{s-1}& \mathrm{if}{\overline{N}}_s^p>N_{high}^p\mathrm{where}{k}_p^1\in [0.5,1);\\ k_p^2{w}_p^{s-1}& \mathrm{if}{\overline{N}}_s^p<N_{low}^p\mathrm{where}{k}_p^2\in (1,2];\\ w_p^{s-1}& \mathrm{if}{N}_{low}^p\le {\overline{N}}_s^p\le N_{high}^p,\end{array}$(A8) ).

6. Use ${\hat{I}}_{target}^m(s,{\eta }_m)=[{\hat{N}}_{inf}^m,{\hat{N}}_{sup}^m]$ for backtracking scheme in Algorithm 2.

In the numerical experiments presented in this paper, we have chosen to use the following line search criteria. In the first numerical experiment, we have chosen $I_{target}^m=[0,26]$ and in the second $I_{target}^m=[0,117]$. Note that these numbers for $N_{sup}^m$ are not arbitrary; they are chosen as a percentage of the total number of interior cells. Furthermore, we choose ${\eta }_m=1/5$, S = 200 and $[N_{low}^m,N_{high}^m]=[\frac{1}{4}{N}_{sup},\frac{3}{4}{N}_{sup}]$ for both experiments. Note that here these intervals are constant for each level set function, but other choices are also possible. These are the parameters that define the dynamic intervals described in (A9(A9) $\begin{array}{rl}{\hat{N}}_{low}^m(s,{\eta }_m)& =-\frac{1}{S}{N}_{low}^m{\eta }_{m}s+N_{low}^m,\end{array}$(A9) ), (A10(A10) $\begin{array}{rl}{\hat{N}}_{high}^m(s,{\eta }_m)& =\frac{1}{S}(N_{low}^m-N_{high}^m){\eta }_{m}s+N_{high}^m;\end{array}$(A10) ).