# Embedding Approach to Modeling Electromagnetic Fields in a Complex Two-Dimensional Environment

Authors: Anton Tijhuis; Ann Franchois; Jean-Michel Geffrin
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

An approach is presented to combine the response of a two-dimensionally inhomogeneous dielectric object in a homogeneous environment with that of an empty inhomogeneous environment. This allows an efficient computation of the scattering behavior of the dielectric cylinder with the aid of the CGFFT method and a dedicated extrapolation procedure. Since a circular observation contour is adopted, an angular spectral representation can be employed for the embedding. Implementation details are discussed for the case of a closed 434 MHz microwave scanner, and the accuracy and efficiency of all steps in the numerical procedure are investigated. Guidelines are proposed for choosing computational parameters such as truncation limits and tolerances. We show that the embedding approach does not increase the CPU time with respect to the forward problem solution in a homogeneous environment, if only the fields on the observation contour are computed, and that it leads to a relatively small increase when the fields on the mesh are computed as well.


---

## Body

## 1. Introduction

In almost any computational approach to solving nonlinear inverse-scattering problems, a discretized configuration is introduced that depends on a fixed number of parameters. Subsequently, a cost functional is defined in terms of simulated and known scattered fields. Here, two different strategies can be distinguished. Conventionally, the corresponding forward problem is treated as an auxiliary problem, which is solved exactly for successive approximate configurations [1–6]. For multidimensional problems, this requires a number of field computations for a varying physical parameter such as frequency or source position. The cost function then refers to the known measured field information and preferably includes a regularizing function of the configuration parameters [7–9]. In the so-called modified gradient method and subsequent generalizations, the configuration and the unknown fields are determined simultaneously [10, 11]. The conventional approach has the advantage that the formulation of the inverse problem directly relates the parameterized configuration to the known field data. From a practical point of view, however, it is sometimes considered as less feasible because of the computational effort required in the repeated field computations. The argument is that it is not needed to compute the field with full accuracy in a configuration that still deviates considerably from the actual one.For the case of an inhomogeneous, lossy dielectric cylinder in a homogeneous surrounding medium, however, it was demonstrated that a highly efficient implementation is obtained when the fields are computed by solving a contrast-source integral equation with a combination of the conjugate-gradient FFT (CGFFT) method and a special extrapolation procedure [12]. The extrapolation can be performed for almost any physical parameter [13], such as frequency or source position. Thus, the forward scattering problem can be solved for each new value of the physical parameter in a few iterations of the CGFFT procedure. This technique has been demonstrated successfully in the context of Newton-type inverse scattering [5, 6]. It is the authors’ experience that these schemes, with only the parameters in the profile parameterization as fundamental unknowns, are generally more efficient than schemes where the field and the profile are determined simultaneously.A special feature of our implementation is that its efficiency is based on the circumstance that the dielectric cylinder is embedded in a homogeneous surrounding medium. This means that Green’s function in the integral equation above exhibits convolution symmetry. Preserving that symmetry in the relevant space discretization allows the application of FFT operations in evaluating the operator products in the conjugate gradient method [14–16]. In practical experiments, however, the surrounding medium may be inhomogeneous and the symmetry is broken. In that case, FFT operations are no longer applicable. The same problem also arises in the modified gradient method, where FFT operations are used as well to compute the field updates.To circumvent this problem, we use the feature that the scattering operator characterizes the complete electromagnetic response of the region inside a closed observation contour. Hence, it must be possible to determine the scattered-field data from a cylinder in an arbitrary environment from the scattering operator for the same object in a homogeneous environment. That data, in turn, can then be obtained with the existing implementation. We introduced this so-called embedding approach for a cylinder inside a circular observation contour in [17, 18]. The choice of this particular configuration was inspired by the experimental research with a 434 MHz scanner of the third author [19]. Several authors have shown interest in quantitative imaging with a circular scanner with metallic enclosure, including [20–22]. It is well-known that employing the 2D Green’s function of the empty casing is computationally expensive [19, 20].In the present paper, we formulate the embedding approach in the angular spectral representation for a general surrounding 2D medium and subsequently specialize to the case where the dielectric cylinder is surrounded by a perfectly conducting circular container. Besides a more comprehensive theoretical formulation than in [17, 18], we provide details on the numerical implementation and on its performance in accuracy and speed as a function of various parameters. We show that, with well-chosen values for these parameters, the embedding approach does not increase the CPU time as compared to the forward problem solution in a homogeneous environment, if only the fields on the observation contour are computed, and that it leads to a relatively small increase in CPU time, when the fields in the object are needed as well, for example, to compute the Jacobian matrix in a Newton-type inversion scheme.The embedding approach relies on the identification of the scattering and reflection operators for the dielectric cylinder and the empty microwave scanner, respectively. The idea of using such an operator to characterize scattering properties has a long tradition in the electromagnetic literature [23–26]. The use of a numerically computed scattering operator, however, is new and originates from the availability of the “march in source position” method [13]. A generalization for multiple interacting domains of arbitrary shapes is given in [27, 28].Finally, it should be remarked that in [29], a procedure is proposed based on reciprocity that is capable of converting the field in the complete configuration into the field in a homogeneous environment, that is, the reverse procedure from what is proposed in the present paper. The suggestion is to perform the profile inversion on the thus corrected data, using inverse-profiling algorithms for objects in a homogeneous background. However, this idea has two possible drawbacks. First, in order to carry out the conversion from one environment to another, complete data on a contour surrounding the scatterer must be available. In an actual experiment, such data may not always be available, while theoretical results for an estimated configuration can always be computed. Second, the conversion renormalizes the experimental data including the measurements, while the present procedure allows a comparison with the actual data. This makes it easier to account for the accuracy of these data, for example, by including appropriate weighting coefficients in the cost functional.The paper is organized as follows. In Section2, we describe the scanner configuration and its mathematical idealization. Section 3 summarizes the field computation for an object in a homogeneous environment. The scattering operators are introduced in Section 4 and used to formulate the embedding approach in Section 5. Section 6 presents the computational details for a homogeneous environment, the empty scanner and the complete configuration. In Section 7, the computational complexity of the algorithm is analyzed, and a procedure is given for tuning the computational parameters. The conclusions are drawn in Section 8.


## 2. Formulation of the Problem

In the present paper, we describe and investigate an efficient procedure to calculate the electric field inside a cylindrical scanner for a given permittivity profile. This field may then be used in inverse-profiling algorithms. For our numerical experiments, we adopted the configuration of the scanner described in [19] and shown in Figure 1, which was developed to conduct biomedical imaging experiments. This scanner comprises a circular array of 64 transmitting/receiving conical dipole antennas, which operate at 434 MHz in a multi-incidence mode, that is, one antenna at a time is transmitting and the others are receiving. The array has a radius of 27.6 cm and is placed inside a water-filled metal casing with a slightly larger radius of 29.0 cm. Measurements of the relative permittivity of the water typically yielded ε1r=76.3+i3.9, which corresponds to a wavelength λ1=7.9 cm in the water. The diameter of the T/R circle thus is about 7λ1, the antennas are spaced apart about λ1/3, and they are at a distance of about λ1/4 from the casing.Figure 1
Synthetic dielectric object in the scanner.Our modeling assumptions are that we may assume the fields to be two-dimensional and that mutual coupling between the antenna elements can be neglected. We therefore consider the configuration shown in Figure2. An inhomogeneous, lossy dielectric cylinder in an observation domain DO is excited by a time-harmonic electric line source on a circular contour ∂DO with radius ρO outside the cylinder. “Complete” scattering data need to be determined, that is, the scattered electric field must be calculated on a circle ∂DO for a line source anywhere on that contour. The “environment” in D¯O is two-dimensional and linearly reacting, that is, the reflected field due to sources in DO∪∂DO satisfies the superposition principle. The observation contour is located in a homogeneous region a<ρ<b with relative permittivity ε1r. The radii a and b bound the cylinder and the environment from the outside and inside, respectively. For now, we leave the environment unspecified, but in the calculations, it will be either a homogeneous space or a conducting cylinder.Figure 2
Two-dimensional model of object in scanner.The electrically polarized field caused by a line source at source pointρS can be identified as Green’s function, that is, the solution of the second-order differential equation
(1)∇T2−s2c02εrρGρ,ρS=−δρ−ρSfor ρ<b that satisfies the proper conditions in ρ>b. In (1), ρ is a two-dimensional position vector and s is a complex frequency with Res≥0. In the present paper, we describe an approach for solving (1) for the case where the scattering configuration is located inside a general environment. This solution will be denoted as Gcasρ,ρS. In our approach, we first determine Green’s function for the same configuration in a homogeneous environment, which will be denoted as Ghomρ,ρS. The differential equation (1) applies to both problems. For the object in a homogeneous environment, Ghomρ,ρS satisfies the radiation condition as ρ→∞. For the object in a general environment, the boundary condition depends on the properties of the exterior medium. For the special case where the environment is a metal wall with inner radius ρ=b, we have Gcasρ,ρS=0 at that wall. For the general case, a reflection operator will be introduced further on that can be employed to formulate a boundary condition on the observation contour ∂DO.


## 3. Homogeneous Environment

The feasibility of our approach depends strongly on the availability of a fast procedure for determining the field for the scattering configuration in a homogeneous environment. Such a scheme is available from [12]. We briefly repeat its main features. Starting point is the contrast-source integral equation
(2)Ghomρ,ρS=G1ρ,ρS−s2c02∬Dχρ′G1ρ,ρ′Ghomρ′,ρSdAρ′,where χρ=εrρ−ε1r, D is a finite domain in which χρ≠0, and G1ρ,ρ′ is Green’s function of the surrounding homogeneous dielectric medium:
(3)G1ρ,ρ′=12πK0sc1ρ−ρ′.In (3), K0 denotes the modified Bessel function of the second kind of order zero. In this integral equation, the incident field and Green’s function are available in closed form.The square region−ℓ<x<ℓ and −ℓ<y<ℓ, in which D is enclosed, is subdivided into N×N subregions with mesh size h=2ℓ/N. The grid points of the square mesh are located at ρm,n=xmux+ynuy, with xm=−ℓ+mh for m=0,1,…,N, and yn=−ℓ+nh for n=0,1,…,N. Solving (2) now amounts to determining an approximation of Ghomρ,ρs at the grid points ρ=ρm,n.The space discretization of the integral in the right-hand side of (2) has two special aspects. First, the logarithmically singular behavior of K0sR/c1 as R=ρ−ρ′↓0 is subtracted by breaking up the integral over D into
(4)∬DK0sRc1+lnRaχρ′Ghomρ′,ρSdAρ′−∬DlnRaχρ′Ghomρ′,ρSdAρ′.Second, the discretization of the integrals in (4) is based on approximating suitable parts of both integrands by piecewise-linear approximations and integrating analytically over polygons determined by the boundary of D and the grid. This results in a discretized integral equation of the form
(5)Gm,n=G1ρm,n,ρS−s2h22πc02∑m′=0N∑n′=0NKm−m′,n−n′χm′,n′Gm′,n′,where χm′,n′ is a sampled, filtered version of χρ. In (5), the convolution-type structure of the continuous equation (2) has been preserved. This makes this equation suitable for the application of the CGFFT method. In addition, it is second-order accurate in the mesh size h. More information on the discretization and the corresponding error estimate can be found in [12].The initial estimate for the CGFFT procedure is obtained by taking a linear combination of previous “final” solutions and determining the coefficients by minimizing the squared error for the problem at hand. In [13], a more detailed explanation is given, as well as several examples of physical parameters for which the effectiveness of this extrapolation has been demonstrated. In the present context, we extrapolate in source position [6]. Since ρS=ρO, the physical parameter that is varied here is the angle φS, which explains the denomination “marching on in angle.”


## 4. Scattering Operators

Before we consider the complete problem, we first introduce the scattering operators for the individual building blocks of the configuration. Since both the observation contour and the casing have circular symmetry, it is convenient to carry out a Fourier transformation with respect to the angleφ and to carry out the analysis in the spectral domain. This leads to a matrix formulation in terms of the angular spectral coefficients.


### 4.1. Scattering by a Dielectric Cylinder

We start by considering the dielectric cylinder. For a general excitation outsideDO, we can write the field in a<ρ<ρO in spectral form as
(6)Ezρ=∑m=−∞∞AmexpimφImsρc1+∑m=−∞∞BmexpimφKmsρc1.The coefficientsAm represent the incident field and the coefficients Bm the field scattered by the dielectric object. Since the cylinder is linearly reacting, these coefficients must be related by a linear scattering operator:
(7)Bm=∑m′=−∞∞Sm,m′Am′.The value of the elementsSm,m′ in (7) can be obtained from Green’s function for the cylinder in a homogeneous environment. In the special case described in Section 3, the incident field is the field caused by a line source in a completely homogeneous background and can be expressed in spectral form as
(8)Ghomincρ,ρS=12π∑m=−∞∞expimφ−φSImsρ<c1Kmsρ>c1,where ρ<=minρ,ρO and ρ>=maxρ,ρO. This expression is valid for all ρ∈ℝ2 , that is, for 0≤ρ<∞.Now, fora<ρ<ρO, the incident fields specified in (6) and (8) must be identical. Since the expansion functions expimφ are linearly independent, comparing both expressions directly leads to the identification
(9)Am=12πKmsρOc1exp−imφS.Using the definition (7) then leads to the following expression for the scattered field for a receiver position ρR∈∂DO for a cylinder in a homogeneous environment
(10)GhomscatρR,ρS=12π∑m=−∞∞∑m′=−∞∞Sm,m′expimφR−im′φSKmsρOc1Km′sρOc1.This means that the individual elementsSm,m′ that define the scattering operator can be obtained from GhomscatρR,ρS by applying Fourier transformations with respect to φR and φS. The modified Bessel function of the second kind Kmp has no zeros in the right half-plane Rep≥0. Therefore, no problems with division by such functions will be encountered.


### 4.2. Reflection by a Cylindrical Casing

Next, we consider the environment. Writing the field inρO<ρ<b in the spectral form (6), where the coefficients Bm correspond to a field radiated by sources in DO and where the coefficients Am represent a source-free field reflected by the environment, leads to the introduction of the reflection operator:
(11)Am=∑m′=−∞∞Rm,m′Bm′.This operator can be determined from the field that is excited by a line source in a configuration where the actual environment inD¯O surrounds an “empty” observation domain DO that contains a homogeneous dielectric medium with ε=ε1.In the computations, we will consider the special case where the environment is a perfectly conducting circular container with inner radiusρ=b. In that case, it follows from the boundary condition Ez=0 that Rm,m′ is a diagonal operator:
(12)Rm,m′=−Kmsb/c1Imsb/c1δm,m′=Rmδm,m′,where δm,m′ is the Kronecker symbol. The function Imsb/c1 can have zeros on the imaginary s-axis. Such a zero corresponds to the situation where the empty casing has an interior mode. However, interior modes occur only when the medium inside the casing is lossless. In the experimental setup, the water in which the dielectric object is immersed is lossy. In the present study, we have therefore not considered special measures to handle interior modes.


## 4.1. Scattering by a Dielectric Cylinder

We start by considering the dielectric cylinder. For a general excitation outsideDO, we can write the field in a<ρ<ρO in spectral form as
(6)Ezρ=∑m=−∞∞AmexpimφImsρc1+∑m=−∞∞BmexpimφKmsρc1.The coefficientsAm represent the incident field and the coefficients Bm the field scattered by the dielectric object. Since the cylinder is linearly reacting, these coefficients must be related by a linear scattering operator:
(7)Bm=∑m′=−∞∞Sm,m′Am′.The value of the elementsSm,m′ in (7) can be obtained from Green’s function for the cylinder in a homogeneous environment. In the special case described in Section 3, the incident field is the field caused by a line source in a completely homogeneous background and can be expressed in spectral form as
(8)Ghomincρ,ρS=12π∑m=−∞∞expimφ−φSImsρ<c1Kmsρ>c1,where ρ<=minρ,ρO and ρ>=maxρ,ρO. This expression is valid for all ρ∈ℝ2 , that is, for 0≤ρ<∞.Now, fora<ρ<ρO, the incident fields specified in (6) and (8) must be identical. Since the expansion functions expimφ are linearly independent, comparing both expressions directly leads to the identification
(9)Am=12πKmsρOc1exp−imφS.Using the definition (7) then leads to the following expression for the scattered field for a receiver position ρR∈∂DO for a cylinder in a homogeneous environment
(10)GhomscatρR,ρS=12π∑m=−∞∞∑m′=−∞∞Sm,m′expimφR−im′φSKmsρOc1Km′sρOc1.This means that the individual elementsSm,m′ that define the scattering operator can be obtained from GhomscatρR,ρS by applying Fourier transformations with respect to φR and φS. The modified Bessel function of the second kind Kmp has no zeros in the right half-plane Rep≥0. Therefore, no problems with division by such functions will be encountered.


## 4.2. Reflection by a Cylindrical Casing

Next, we consider the environment. Writing the field inρO<ρ<b in the spectral form (6), where the coefficients Bm correspond to a field radiated by sources in DO and where the coefficients Am represent a source-free field reflected by the environment, leads to the introduction of the reflection operator:
(11)Am=∑m′=−∞∞Rm,m′Bm′.This operator can be determined from the field that is excited by a line source in a configuration where the actual environment inD¯O surrounds an “empty” observation domain DO that contains a homogeneous dielectric medium with ε=ε1.In the computations, we will consider the special case where the environment is a perfectly conducting circular container with inner radiusρ=b. In that case, it follows from the boundary condition Ez=0 that Rm,m′ is a diagonal operator:
(12)Rm,m′=−Kmsb/c1Imsb/c1δm,m′=Rmδm,m′,where δm,m′ is the Kronecker symbol. The function Imsb/c1 can have zeros on the imaginary s-axis. Such a zero corresponds to the situation where the empty casing has an interior mode. However, interior modes occur only when the medium inside the casing is lossless. In the experimental setup, the water in which the dielectric object is immersed is lossy. In the present study, we have therefore not considered special measures to handle interior modes.


## 5. The Complete Configuration

Now that the individual building blocks of the configuration have been analyzed, we can combine the results to obtain the field in the complete configuration.


### 5.1. Field on the Observation Contour

We first determine the field on the observation contour∂DO. To that end, we divide the domain a<ρ<b, in which the spectral analysis is carried out, in two subdomains that are separated by ∂DO.In the complete regiona<ρ<b, we have the spectral representation
(13)Gcasρ,ρS=12π∑m=−∞∞expimφ−φSImsρ<c1Kmsρ>c1+∑m=−∞∞expimφCmImsρc1+∑m=−∞∞expimφDmKmsρc1,where the coefficients Cm and Dm depend on the angle φS.For the “interior subdomain”a<ρ≤ρO, we may envisage the situation as scattering by the dielectric cylinder. We treat the field due to the line source and the field reflected from the casing (given by the coefficients Cm) as primary fields that generate a secondary field (given by the coefficients Dm) that propagates in the direction of increasing ρ. With the aid of the definition (7), we then directly arrive at
(14)Dm=∑m′=−∞∞Sm,m′exp−im′φS2πKm′sρOc1+Cm′.For the “exterior subdomain”ρO≤ρ<b, we similarly treat the line source and the field scattered by the dielectric cylinder as primary fields that generate a secondary field reflected from the metal casing. This leads to
(15)Cm=∑m′=−∞∞Rm,m′exp−im′φS2πIm′sρOc1+Dm′.Combining both results finally gives a linear equation for the coefficientsCm:
(16)Cm−∑m′=−∞∞Rm,m′∑m″=−∞∞Sm′,m″Cm″=12π∑m′=−∞∞Rm,m′exp−im′φSIm′sρOc1+12π∑m′=−∞∞Rm,m′∑m″=−∞∞Sm′,m″exp−im″φSKm″sρOc1.This equation can be solved by truncating the summation overm′ and m″ and inverting the resulting matrix equation. This must be carried out for varying φS, that is, for multiple right-hand sides. Once the coefficients Cm are found, we can use expression (14) to determine Dm. Substitution of the values for these coefficients in (13) then gives the total field Gcasρ,ρS in a<ρ<b.


### 5.2. Operator Formulation

The equation found in (16) can also be written in operator form. From a spectral point of view, the line source generates a source-free incident field in ρ<ρO and a radiating field in ρ>ρO, with spectral amplitudes
(17)Am<=12πKmsρOc1exp−imφS,Am>=12πImsρOc1exp−imφS,respectively. Identifying these amplitudes as the elements of two excitation vectors A< and A> then results in the operator equation
(18)1−RSC=RA>+RSA<,where the unknown vector C contains the spectral amplitudes for the secondary field that originates from the casing. From energy considerations, it follows that the norm of the operator product RS is sufficiently small that the solution of (18) may formally be written as a geometrical series:
(19)C=∑n=0∞RSnRA>+RSA<.This confirms that repeated scattering and diffraction effects are accounted for in solving the system of equations given by (16) and (18). Moreover, (19) constitutes the basis for a Neumann-type iterative scheme that can be used for the solution of these equations. In the same notation, we have from (14)
(20)D=SA<+C,which completes the formulation in operator form.


### 5.3. Field in the Observation Domain

The previous analysis allows us to compute the fieldGcasρR,ρS for ρR and ρS on ∂DO. However, in inverse-scattering algorithms, we must also determine Gcasρ,ρS for ρ∈DO and for ρS∈∂DO [5, 18]. This field is needed to determine the “profile update” in Newton-type optimization.To this end, we use the equivalence principle. InDO, the total field may be envisioned as a response to the current at the line source and the induced surface current on the casing. Both these currents radiate a field that is incident on the dielectric cylinder. Alternatively, we may treat the second constituent as originating from an equivalent surface current on ∂DO. This means that we replace the right-hand side of (1) by
(21)−δρ−ρS−∫02πwφP,φSδρ−ρPdφPand compute the field generated by this source for a dielectric cylinder in a homogeneous environment. In DO, we will then find the correct field Gcasρ,ρS. This conclusion is in fact merely a special formulation of Huygens’ principle. Obviously, the field in domain D¯O of this equivalent configuration will deviate from the actual field.Determining the incident field by separation of variables and comparing the result fora<ρ<ρO with the terms involving Imsρ/c1 in (13) lead to the definition
(22)wφP,φS=∑m=−∞∞CmφSKmsρO/c1expimφP.Now, the uniqueness of the solution of the homogeneous wave equation implies that, for a given incident field inDO, the corresponding total field in DO is fixed. Therefore, we may identify the total field in DO as a superposition of fields generated by line sources on ∂DO in the homogeneous embedding
(23)Gcasρ,ρS=Ghomρ,ρS+∫02πwφP,φSGhomρ,ρPdφP,which of course holds only for ρ<ρO. In (23), ρP is a point on ∂DO, characterized by the angle φP. With (23), we can now compute Gcasρ,ρS directly from Ghomρ,ρP with 0<φP<2π.For a general environment, an analogous procedure is available for the field inD¯O. We then express the field in the environment in terms of the line-source response of the empty casing. For a circular metal casing, such a procedure is not needed, since the spectral representation (13) remains valid up to ρ=b.


### 5.4. Direct Reflection from the Casing

In the practical configuration shown in Figure1, the line source is located close to the interface at ρ=b. This means that the direct response, corresponding to the term RA> in (19), must compensate the logarithmically singular behavior of the incident field as ρ−ρS↓0. In fact, for φ≈φS and ρ≈b, the field directly reflected from the casing may be approximated by that of an image source with opposite sign at ρ≈2b−ρO and φ=φS. Therefore, the convergence of the angular series slows down as ρO↑b. Now, it is observed from (12) that the computation of the reflection coefficients Rm is much easier than the computation of the elements of the scattering matrix Sm,m′. Therefore, it makes sense to extract the direct reflection from the casing, which gives rise to the term RA>, out of the total casing response determined by the coefficients Cm. To this end, we rewrite the operator equation as
(24)1−RSC−RA>=RSRA>+A<,which leads to the power-series solution
(25)C~=C−RA>=∑n=1∞RSnRA>+A<.The solution (25) separates the field that originates from the casing into a field that would be present in an empty casing and an additional field scattered by the dielectric. Both the field excited by the line source and the field that results from a direct reflection by the casing are considered as fields that are incident on the dielectric cylinder. This combined incident field is scattered at least once by the dielectric cylinder and reflected at least once at the casing before it contributes to the regularized spectral coefficient C~. Since the distance between this cylinder and the observation contour is considerably larger than the distance between the casing and the observation contour, the resulting diffracted field is much smoother than the fields that originate from the line source and from direct reflection at the casing.To facilitate the discussion of the numerical aspects in the upcoming sections, we express the operator form of the decomposition (25) of the field in the region a<ρ<b in the spectral form used in Section 5.1. Further, we restrict ourselves to the diagonal operator specified in (12). Combining (17) and (24) leads to the form
(26)Gcasρ,ρS=12π∑m=−∞∞expimφ−φSImsρ<c1Kmsρ>c1+∑m=−∞∞RmAm>expimφImsρc1+∑m=−∞∞C~mexpimφImsρc1+∑m=−∞∞DmexpimφKmsρc1.The second sum in (26) represents the field that is directly reflected from the casing, with Rm given by (12) and Am> given by (17). The field Gcasρ,ρS can also be written as the sum
(27)Gcasρ,ρS=Gcasemptyρ,ρS+Gcasdifρ,ρS,where Gcasemptyρ,ρS is the field in the empty casing—the first two sums in (26)—and Gcasdifρ,ρS is the difference field due to the insertion of the object—the last two sums in (26). In an experiment, GcasdifρR,ρS is obtained by subtracting the measured fields with and without the object in place.In a similar fashion, the equivalent surface current in (21) is decomposed as
(28)wφP,φS=wemptyφP,φS+wdifφP,φS,where
(29)wemptyφP,φS=12π∑m=−∞∞expimφ−φSRImmsρO/c1KmsρO/c1corresponds to the direct reflection from the casing and where wdifφP,φS is obtained by replacing CmφS by C~mφS in (22). wdifφP,φS generates the source-free constituent of the difference field, that is, the third sum in (26).


## 5.1. Field on the Observation Contour

We first determine the field on the observation contour∂DO. To that end, we divide the domain a<ρ<b, in which the spectral analysis is carried out, in two subdomains that are separated by ∂DO.In the complete regiona<ρ<b, we have the spectral representation
(13)Gcasρ,ρS=12π∑m=−∞∞expimφ−φSImsρ<c1Kmsρ>c1+∑m=−∞∞expimφCmImsρc1+∑m=−∞∞expimφDmKmsρc1,where the coefficients Cm and Dm depend on the angle φS.For the “interior subdomain”a<ρ≤ρO, we may envisage the situation as scattering by the dielectric cylinder. We treat the field due to the line source and the field reflected from the casing (given by the coefficients Cm) as primary fields that generate a secondary field (given by the coefficients Dm) that propagates in the direction of increasing ρ. With the aid of the definition (7), we then directly arrive at
(14)Dm=∑m′=−∞∞Sm,m′exp−im′φS2πKm′sρOc1+Cm′.For the “exterior subdomain”ρO≤ρ<b, we similarly treat the line source and the field scattered by the dielectric cylinder as primary fields that generate a secondary field reflected from the metal casing. This leads to
(15)Cm=∑m′=−∞∞Rm,m′exp−im′φS2πIm′sρOc1+Dm′.Combining both results finally gives a linear equation for the coefficientsCm:
(16)Cm−∑m′=−∞∞Rm,m′∑m″=−∞∞Sm′,m″Cm″=12π∑m′=−∞∞Rm,m′exp−im′φSIm′sρOc1+12π∑m′=−∞∞Rm,m′∑m″=−∞∞Sm′,m″exp−im″φSKm″sρOc1.This equation can be solved by truncating the summation overm′ and m″ and inverting the resulting matrix equation. This must be carried out for varying φS, that is, for multiple right-hand sides. Once the coefficients Cm are found, we can use expression (14) to determine Dm. Substitution of the values for these coefficients in (13) then gives the total field Gcasρ,ρS in a<ρ<b.


## 5.2. Operator Formulation

The equation found in (16) can also be written in operator form. From a spectral point of view, the line source generates a source-free incident field in ρ<ρO and a radiating field in ρ>ρO, with spectral amplitudes
(17)Am<=12πKmsρOc1exp−imφS,Am>=12πImsρOc1exp−imφS,respectively. Identifying these amplitudes as the elements of two excitation vectors A< and A> then results in the operator equation
(18)1−RSC=RA>+RSA<,where the unknown vector C contains the spectral amplitudes for the secondary field that originates from the casing. From energy considerations, it follows that the norm of the operator product RS is sufficiently small that the solution of (18) may formally be written as a geometrical series:
(19)C=∑n=0∞RSnRA>+RSA<.This confirms that repeated scattering and diffraction effects are accounted for in solving the system of equations given by (16) and (18). Moreover, (19) constitutes the basis for a Neumann-type iterative scheme that can be used for the solution of these equations. In the same notation, we have from (14)
(20)D=SA<+C,which completes the formulation in operator form.


## 5.3. Field in the Observation Domain

The previous analysis allows us to compute the fieldGcasρR,ρS for ρR and ρS on ∂DO. However, in inverse-scattering algorithms, we must also determine Gcasρ,ρS for ρ∈DO and for ρS∈∂DO [5, 18]. This field is needed to determine the “profile update” in Newton-type optimization.To this end, we use the equivalence principle. InDO, the total field may be envisioned as a response to the current at the line source and the induced surface current on the casing. Both these currents radiate a field that is incident on the dielectric cylinder. Alternatively, we may treat the second constituent as originating from an equivalent surface current on ∂DO. This means that we replace the right-hand side of (1) by
(21)−δρ−ρS−∫02πwφP,φSδρ−ρPdφPand compute the field generated by this source for a dielectric cylinder in a homogeneous environment. In DO, we will then find the correct field Gcasρ,ρS. This conclusion is in fact merely a special formulation of Huygens’ principle. Obviously, the field in domain D¯O of this equivalent configuration will deviate from the actual field.Determining the incident field by separation of variables and comparing the result fora<ρ<ρO with the terms involving Imsρ/c1 in (13) lead to the definition
(22)wφP,φS=∑m=−∞∞CmφSKmsρO/c1expimφP.Now, the uniqueness of the solution of the homogeneous wave equation implies that, for a given incident field inDO, the corresponding total field in DO is fixed. Therefore, we may identify the total field in DO as a superposition of fields generated by line sources on ∂DO in the homogeneous embedding
(23)Gcasρ,ρS=Ghomρ,ρS+∫02πwφP,φSGhomρ,ρPdφP,which of course holds only for ρ<ρO. In (23), ρP is a point on ∂DO, characterized by the angle φP. With (23), we can now compute Gcasρ,ρS directly from Ghomρ,ρP with 0<φP<2π.For a general environment, an analogous procedure is available for the field inD¯O. We then express the field in the environment in terms of the line-source response of the empty casing. For a circular metal casing, such a procedure is not needed, since the spectral representation (13) remains valid up to ρ=b.


## 5.4. Direct Reflection from the Casing

In the practical configuration shown in Figure1, the line source is located close to the interface at ρ=b. This means that the direct response, corresponding to the term RA> in (19), must compensate the logarithmically singular behavior of the incident field as ρ−ρS↓0. In fact, for φ≈φS and ρ≈b, the field directly reflected from the casing may be approximated by that of an image source with opposite sign at ρ≈2b−ρO and φ=φS. Therefore, the convergence of the angular series slows down as ρO↑b. Now, it is observed from (12) that the computation of the reflection coefficients Rm is much easier than the computation of the elements of the scattering matrix Sm,m′. Therefore, it makes sense to extract the direct reflection from the casing, which gives rise to the term RA>, out of the total casing response determined by the coefficients Cm. To this end, we rewrite the operator equation as
(24)1−RSC−RA>=RSRA>+A<,which leads to the power-series solution
(25)C~=C−RA>=∑n=1∞RSnRA>+A<.The solution (25) separates the field that originates from the casing into a field that would be present in an empty casing and an additional field scattered by the dielectric. Both the field excited by the line source and the field that results from a direct reflection by the casing are considered as fields that are incident on the dielectric cylinder. This combined incident field is scattered at least once by the dielectric cylinder and reflected at least once at the casing before it contributes to the regularized spectral coefficient C~. Since the distance between this cylinder and the observation contour is considerably larger than the distance between the casing and the observation contour, the resulting diffracted field is much smoother than the fields that originate from the line source and from direct reflection at the casing.To facilitate the discussion of the numerical aspects in the upcoming sections, we express the operator form of the decomposition (25) of the field in the region a<ρ<b in the spectral form used in Section 5.1. Further, we restrict ourselves to the diagonal operator specified in (12). Combining (17) and (24) leads to the form
(26)Gcasρ,ρS=12π∑m=−∞∞expimφ−φSImsρ<c1Kmsρ>c1+∑m=−∞∞RmAm>expimφImsρc1+∑m=−∞∞C~mexpimφImsρc1+∑m=−∞∞DmexpimφKmsρc1.The second sum in (26) represents the field that is directly reflected from the casing, with Rm given by (12) and Am> given by (17). The field Gcasρ,ρS can also be written as the sum
(27)Gcasρ,ρS=Gcasemptyρ,ρS+Gcasdifρ,ρS,where Gcasemptyρ,ρS is the field in the empty casing—the first two sums in (26)—and Gcasdifρ,ρS is the difference field due to the insertion of the object—the last two sums in (26). In an experiment, GcasdifρR,ρS is obtained by subtracting the measured fields with and without the object in place.In a similar fashion, the equivalent surface current in (21) is decomposed as
(28)wφP,φS=wemptyφP,φS+wdifφP,φS,where
(29)wemptyφP,φS=12π∑m=−∞∞expimφ−φSRImmsρO/c1KmsρO/c1corresponds to the direct reflection from the casing and where wdifφP,φS is obtained by replacing CmφS by C~mφS in (22). wdifφP,φS generates the source-free constituent of the difference field, that is, the third sum in (26).


## 6. Numerical Study

We discuss the numerical implementation of the embedding approach, and we examine its performance in accuracy and speed as a function of various parameters for the configuration of the 434 MHz scanner [19]. All simulations were performed on a SUN ULTRA HPC 4000 workstation. In our code, we have used the public domain software packages LAPACK [30], for the linear system solutions; AMOS [31], for the Bessel function computations; and NMS [32], for the 2D FFTs. In particular, the accuracy of the various building blocks in the embedding approach is evaluated with tests on a number of homogeneous circular dielectric cylinders, specified in Table 1, for which accurate field solutions are available from spectral representations. This accuracy is quantified by means of the normalized root-mean-square error (NRMSE), defined as
(30)NRMSE=∑juj−ujexact2∑jujexact2.Table 1
Diameter, relative permittivity and mesh size of the test objects: homogeneous dielectric circular cylinders Muscle1, Muscle4, Air1 and Air2, and an inhomogeneous dielectric cylinder Leg2. The permittivity for muscle is from [19], that for bone and fat is based on [33].
(a)
Cylinder
Muscle1
Muscle4
Air1
Air2
Diameter (cm)
8.8
35.2
8.8
17.6
εr
54.2 + i38.4
54.2 + i38.4
1
1
2ℓ (cm)
9
36
9
18
(b)
Cylinder
Leg2: 3 layers
Fat
Muscle
Bone
Diameter (cm)
17.6
15.6
4.4
εr
6 + i2.5
54.2 + i38.4
6 + i2.5
2ℓ (cm)
18Besides the homogeneous cylinders of muscle and air, with approximate diametersλ1,2λ1, and 4λ1, there is also an inhomogeneous leg object with approximate diameter 2λ1, which consists of a circular cylinder of muscle, covered with a 1 cm thick layer of fat and containing a bone with approximate diameter λ1/2. The bone is decentered 1.5 cm in the x-direction. The homogeneous cylinders are centered at the origin (i.e., the center of the casing), while the leg is decentered by −0.5 cm in the x-direction. On the one hand, the discretization cell size h=2ℓ/N and hence the number of grid points Ngrid=N+1×N+1 determine the accuracy of Ghomρ,ρS on the grid, of GhomscatρR,ρS on ∂DO and of all quantities derived from these in the embedding approach. On the other hand, a sufficient number of line sources—assumed to be equally spaced on ∂DO—are needed to obtain accurate representations for the scattering matrix (10), say M, and for the equivalent surface current on ∂DO (21), say L. The number of forward problem solutions K≥maxL,M may differ from the actual number of sources used in the cost function for an imaging experiment. In this section, we illustrate the influence of K,L,M, and N on the computational performance.


### 6.1. Homogeneous Environment

#### 6.1.1. Discretization Error—Choice ofN

We compared the total fieldGhomρ,ρS in the grid points, computed by solving (5) with the CGFFT method for the source position φS=0, with the exact solution as a function of the discretization cell size h. For the nonsymmetrical Leg2 object, the exact solution was chosen as the discretized solution for a much smaller cell size. We also compared the corresponding scattered field GhomscatρR,ρS in the 64 receiver points with the exact solution. Figure 3(a) shows the discretization error on the total field, with h ranging from 0.0352 cm≈λ1/224 to 2.25 cm ≈λ1/3.5. Depending on the diameter of the object, we used a mesh with side 2ℓ= 9 cm, 18 cm, or 36 cm, and we varied N between 4 and 512 (or Ngrid between 25 and 263,169). Figure 3(b) shows the resulting discretization error on the scattered field. It is clearly illustrated that the discretization errors are of Oh2 as h↓0, and it furthermore appears that for fixed h, the variation in the errors for the different cylinders considered is limited. Let us add here that for highly contrasting objects, such as air cylinders, extremely small cells may be needed in areas where the field rapidly varies—this occurs, for example, when the air interface is in the neighborhood of the source. For objects with characteristics within the ranges considered in Table 1, Figure 3(b) is helpful for choosing the cell size corresponding to a given acceptable error on the scattered field, and Figure 3(a) then gives an indication for the CGFFT stop criterion. Proceeding with the iterations when the CGFFT NRMSE (with respect to the incident field on the mesh) is smaller than a certain fraction of the total field discretization error is a waste of computational effort.Figure 3
The NRMSE for (a)Ghomρ,ρS on the mesh and for (b) GhomscatρR,ρS on ∂DO as a function of the cell size h for the objects of Table 1. In presence of the casing and for proper choices of M and L, the values in (a) are multiplied with a factor between 1 and 2 for the NRMSE of Gcasρ,ρS, and (b) remains valid for the NRMSE of GcasdifρR,ρS.
(a)
(b)


#### 6.1.2. CPU Time

The CPU times for the different steps in the forward problem are illustrated in Figure4 for Leg2 as a function of N. Most of the time is taken up by the CGFFT iterations to solve (5) for the different angles of incidence; only a few percent of this time is needed to compute, respectively, the initial estimates for the total field by means of the “marching on in angle” procedure, the Green functions K in (5), and the K scattered fields GhomscatρR,ρS in the K receivers. The second most time-consuming step is the computation of the K incident fields G1ρ,ρS on the grid; the values presented here correspond to a worst-case scenario, since we did not exploit the symmetries in the grid and T/R points. Note that the Green functions can be stored for a given configuration of the mesh, transmitters/receivers, and exterior medium permittivity.Figure 4
CPU times for the different steps in the forward problem for Leg2 withK=64 sources: (1) solution of equation (5) for Ghomρ,ρS with CGFFT and “marching on in angle”; (2) computation of the incident field G1ρ,ρS; (3) computation of the scattered field GhomscatρR,ρS; (4) one CGFFT iteration; (5) the summation (42) to obtain the total field Gcasρ,ρS. The curve (6) gives the CPU time for step (1) with the incident field as initial guess.


#### 6.1.3. Marching on in Angle

We tested the efficiency of “marching on in angle,” where we used the total field solutions from the three previous excitations to compute the initial estimate, by comparing the total number of CGFFT iterations to those in a conventional CGFFT forward problem solution, where the initial estimate is chosen equal to the incident field. Such a study has not yet been reported for the case of a relatively large value of the exterior medium permittivityε1r. When the contrast εrρ−ε1r/ε1r is small, only a limited number of CGFFT iterations is needed to obtain convergence.In all our tests, “marching on in angle” was faster, except for the case where the spacing between the sources was larger thanλ1. With proper choices for the number of sources K, for N and for the CGFFT stop criterium, “marching on in angle” leads to a significant reduction in the computational effort. For example, with N corresponding to a scattered field discretization error of a few percent (see Figure 3(b)), the relative reduction typically is 20% for Muscle1 N=16,K=32, 52% for Leg2 N=32,K=64—see also Figure 4 where the squares are below the dotted curve—and 32% for Muscle4 N=64,K=128, see Table 2, which gives the relative reduction in the total number of CGFFT iterations as a function of N and K. The gain in efficiency is not as spectacular as with some examples in [13], for which a much larger number of CGFFT iterations is needed to obtain convergence, which leaves more room for improving the computational efficiency.Table 2
Comparison of the total number of CGFFT iterations as a function ofN and K for the conventional (i.e., incident field initial guess) and “marching on in angle” approaches. The relative reduction in the number of iterations, the CGFFT stop criterium (NMRS CGFFT), and the resulting NRMSE for Ghom on the grid also are indicated.
N
16
32
64
128
256
512
Muscle1
K = 16
No marching
92
108
112
136
Marching
96
111
127
139
Reduction (%)
−4
—
—
—
K = 32
No marching
188
220
224
280
312
344
Marching
151
163
216
256
267
320
Reduction (%)
20
26
4
9
14
7
K = 64
No marching
380
444
448
568
632
696
Marching
176
200
325
390
460
518
Reduction (%)
54
55
27
31
27
26
NRMSE CGFFT
2 × 10−3
10−3
2 × 10−4
5 × 10−5
10−5
2 × 10−6
NRMSEGhom
5.9 × 10−3
1.7 × 10−3
4.4 × 10−4
1.1 × 10−4
2.8 × 10−5
7.1 × 10−6
Leg2
K = 32
No marching
1043
1272
1571
Marching
780
1015
1406
Reduction (%)
25
20
11
K = 64
No marching
2095
2555
3170
3489
3755
3981
Marching
878
1216
1896
2673
3239
3607
Reduction (%)
58
52
40
23
14
9
K = 128
No marching
4199
5112
6352
7533
Marching
1856
2466
3135
4596
Reduction (%)
56
52
49
39
NRMSE CGFFT
10−2
5 × 10−3
10−3
2 × 10−4
5 × 10−5
10−5
NRMSEGhom
1.1 × 10−1
3.2 × 10−2
7.9 × 10−3
1.9 × 10−3
5.1 × 10−4
1.2 × 10−4
Muscle4
K = 64
No marching
700
832
1128
1376
1648
1752
Marching
619
761
1032
1302
1563
1678
Reduction (%)
12
9
9
10
5
4
K = 128
No marching
1404
1664
2256
2760
3312
Marching
791
911
1531
2000
2586
Reduction (%)
44
45
32
28
22
K = 256
No marching
2804
3328
4512
5536
Marching
1034
1103
1679
2565
Reduction (%)
63
67
63
54
NRMSE CGFFT
2 × 10−2
10−2
2 × 10−3
5 × 10−4
10−4
5 × 10−5
NRMSEGhom
6.4 × 10−2
2.6 × 10−2
6.7 × 10−3
1.7 × 10−3
4.5 × 10−4
1.3 × 10−4It appears that the relative reduction decreases with increasing values ofN, that is, finer meshing. This can be explained as follows. The dimension of the vector space in which the solution needs to be determined increases as ON2, while the discretization error, hence the recommended CGFFT stop criterion, decreases as ON−2. The error on the “marching on in angle” initial estimate does not depend on N, while the number of iterations, which is needed to further reduce the error, increases and becomes less dependent on the choice of the initial guess.It can furthermore be seen that the relative reduction increases with increasing numbers of sourcesK. For example, when K is doubled in the aforementioned examples, the relative reduction typically is 54% for Muscle1, 52% for Leg2, and 63% for Muscle4. The sources are closer, such that the error in the “marching on in angle” initial estimate is reduced. This is also illustrated in Figure 5, which shows the number of CGFFT iterations, with and without “marching on in angle” as a function of K for Muscle1.Figure 5
The total number of CGFFT iterations as a function of the number of sourcesK for the “marching on in angle” and conventional approaches, for Muscle1 with N=16.Let us conclude by stressing the importance of choosing the largest possible value for the CGFFT stop criterion (NRMSE CGFFT), in order to get the most benefit from “marching on in angle.” When, for example, for Leg2N=32,K=64, the CGFFT stop criterion is reduced by a factor of 10, the number of CGFFT iterations is almost doubled, from 1216 to 2148, while the resulting reduction in the NRMSE on the total field, from 3.2 to 2.8%, is not significant.


#### 6.1.4. Scattering Matrix—Choice ofM

For the computation of the scattering matrix, (10) is expressed in M transmitter/receiver positions on ∂DO, and the summations are truncated to M terms,
(31)GhomscatρR,ρS=12π∑m=−M/2+1M/2∑m′=−M/2+1M/2Sm,m′KmsρOc1Km′sρOc1expi2πmr−m′tM,with r,t=1,…,M, and where M should be large enough for the aliasing effects to be negligible or Sm+lM,m′+l′MKm+lMKm′+l′M≪Sm,m′KmKm′ for ±l,l′=1,…,∞. The elements
(32)Qm,m′=M2πSm,m′KmsρOc1Km′sρOc1are computed from (31) by means of a 2D FFT. The computational effort of this step is negligible: 10 ms for the example of Figure 4. For a centered circular homogeneous cylinder, the scattering matrix is diagonal with Qm,m given by
(33)−M2πc1Insρ2/c1In−1sρ2/c2−c2Insρ2/c2In−1sρ2/c1c1Knsρ2/c1In−1sρ2/c2−c2Insρ2/c2Kn−1sρ2/c1Km2sρOc1,with ρ2 and c2=c0/ϵ2r the radius of and the wave speed in the cylinder, respectively. For high orders m, the elements (33) decrease with increasing order at a rate which primarily depends on ρ2, for a given ρO. This is illustrated in Table 3, which shows that the truncation errors in (31) are negligible when M=32>2ODP for the smallest cylinders, which yields double precision, and M=64>2OSP for the largest cylinder, which is largely sufficient to obtain single precision. The accuracy of the scattering matrix then is determined by the discretization error on Ghomscat.Table 3
The ordersm=OSP and m=ODP for which the elements Qmm of the scattering matrix have dropped to fractions 10−7 (single precision) and 10−14 (double precision), respectively, of their maximum absolute value, for the homogeneous cylinders of Table 1.
Cylinder
Muscle1
Air1
Air2
Muscle4
OSP
9
9
14
24
ODP
14
14
20
36


### 6.2. Direct Reflection by the Casing

The direct reflection from the casing, which is treated as a separate constituent in the expressions for the field on∂DO and for the equivalent surface current, does not depend on the object. Hence, these contributions can be computed once beforehand for a given water-filled scanner geometry. Let us therefore consider the case of a line source on ∂DO in the casing without object. The spectral representation of the field Gcasemptyρ,ρS in the complete region 0<ρ<b is then given by the first two sums in (26).


#### 6.2.1. Exact Solution

The convergence of the first sum in (26) is extremely slow for observation points on ∂DO. This is visible in Figure 6, in Table 4, and from applying the asymptotic approximation ([34], 9.3.1) for high orders m≫∣s/c1ρ>∣,
(34)12πImsρc1KmsρOc1≈14πmρ<ρ>m.Figure 6
Absolute values of the terms in (26) for ρ=ρO: (1) 1/2π∣ImKm∣, (2) 1/2π∣RmImIm∣, ∣C~mIm∣ computed with equation (40) for (3) Muscle1, (4) Leg2, and (5) Muscle4. The terms ∣DmKm∣ are indicated with +.Table 4
The ordersm=OSP and m=ODP for which the terms in the series in (8) and (35) have dropped to fractions 10-7 (single precision) and 10-14 (double precision), respectively, of their maximum absolute values, for different ρ.
ρ (cm)
27.6
25
15
5
(8): ∣Imsρ/c1KmsρO/c1∣
ρ/ρO
1
0.91
0.54
0.18
OSP
≈107
130
30
15
ODP
≈1014
283
52
22
(35): ∣RmImsρ/c1ImsρO/c1∣
ρ/ρ1
0.91
0.82
0.49
0.16
OSP
97
60
27
15
ODP
211
124
44
22Consequently, it is more efficient to computeG1ρ,ρS with (3) or
(35)Gcasemptyρ,ρS=12πK0sc1ρ−ρS+12π∑m=−OOexpimφ−φSRmImsρOc1Imsρc1,where the reflected-field constituent is truncated at order O. The convergence of the spectral representation of the direct casing reflection in (35) is better than that of (8) in the region 0<ρ<ρO; see Figure 6 and Table 4. In this case, the asymptotic approximation of the term m, for m≫∣s/c1b∣, yields
(36)12πRmImsρOc1Imsρc1→−14πmρρOb2m,which is identical to that of a line source located on a contour with radius ρ1=b2/ρO in a homogeneous background—b thus is the geometric mean of ρO and ρ1. We observed that the direct casing reflection reaches stable and accurate DP values in the 64 receivers on ∂DO when choosing O=ODP=211.


#### 6.2.2. Expression from Equivalent Line Sources

Alternatively, in the regionρ<ρO, the field can be regarded as if it originated from an equivalent surface current (21), which we replace with a discrete set of L equivalent line sources with spacing ΔφP=2π/L on ∂DO in a homogeneous background, denoting by WemptyφP,φS the complex amplitude of the equivalent source at φP due to an excitation at φS. Imposing the identity of the reflected field constituent in (35) and the field generated by these sources, for ρ<ρO,
(37)12π∑m=−L/2+1L/2expimφ−φSRmImsρOc1Imsρc1=12π∑φPWemptyφP,φS∑m=−L/2+1L/2expimφ−φp⋅Imsρc1KmsρOc1,where we used the spectral representation (8) for G1ρ,ρP, we obtain the following expression for the complex amplitude
(38)WemptyφP,φS=1L∑m=−L/2+1L/2expimφP−φSRmImsρO/c1KmsρO/c1.According to (23) and since in this section Ghomρ,ρS=G1ρ,ρS, the field in the region ρ<ρO then is expressed as the linear combination
(39)Gcasemptyρ,ρS=G1ρ,ρS+∑φPWemptyφP,φSG1ρ,ρP.The truncation toL terms in the RHS of (37) may lead to significant errors on the field (39) in grid points that are near to ∂DO. In Table 5, we compared the field on the mesh (39), generated by different numbers L of equivalent line sources, with the exact solution (35) for meshes with sides 9 cm, 18 cm, and 36 cm. It follows that for the mesh with side 9 cm, L=32 already yields a very high precision; for the mesh with side 36 cm, L=64 yields a moderate precision, and L≥128 is needed for high precision computations. Figure 7 shows an image of the amplitude and phase of the field (35) in the water-filled casing.Table 5
The NRMSE of the field (39) in the empty casing as a function of the number of equivalent sources L for different sizes 2ℓ of the mesh.
2ℓ (cm)
L=16
L=32
L=64
L=128
L=256
9
5 ×10−4
2 ×10−9
2 ×10−14
2 ×10−14
18
6 ×10−2
7 ×10−5
4 ×10−14
4 ×10−14
36
10−1
2 ×10−4
3 ×10−7
5 ×10−11Figure 7
Amplitude (a) and phase (b) of the field (35) in the empty water-filled casing. An upper limit of 0.2 was chosen in the image of the amplitude.
(a)
(b)


### 6.3. Complete Configuration

Finally, we look into the computation of the fields for the complete configuration of an object in the casing. It is shown that the embedding approach maintains the accuracy of the forward problem solution in homogeneous space, ifM (order of the scattering matrix) and L (number of equivalent line sources) are properly chosen. The additional computational effort for the embedding as such is also examined.


#### 6.3.1. Scattered Field on the Receivers

Instead of solving (16) for Cm, we take into account the separation of the direct casing reflection (26) and solve the resulting set of equations, truncated to M terms, for C~mIm
(40)C~mφSImsρOc1−Um∑n=−M/2+1M/22πVnQmnMC~nφSInsρOc1=Um∑n=−M/2+1M/2QmnM1+Unexp−inφS,where Um=RmIm/Km and Vm=ImKm with m=−M/2+1…M/2. The elements DmKm then are computed with
(41)DmφSKmsρOc1=C~mφSImsρO/c1Um.The absolute values ofC~mIm and DmKm rapidly decrease as a function of m, as is shown in Figure 6. As a consequence, the convergence of the spectral representation (26) of the field GcasdifρR,ρS is much better than that of GcasemptyρR,ρS. Provided that M is chosen according to Table 3, the error on GcasdifρR,ρS only depends on the discretization error in the scattered field GhomscatρR,ρS. We compared the values of GcasdifρR,ρS with the exact solution as a function of N for Muscle1 with M=32, Leg2 with M=64, and Muscle4 with M=128, and we observed that the resulting NRMSE is almost identical to that for the homogeneous solutions; see Figure 3(b).


#### 6.3.2. Total Field on the Grid

The total fieldGcasρ,ρS on the mesh is computed as a linear combination of homogeneous solutions for L equivalent line sources
(42)Gcasρ,ρS=Ghomρ,ρS+∑φPWemptyφP,φS+WdifφP,φSGhomρ,ρP,where WemptyφP,φS is given by (38) and where
(43)WdifφP,φS=2πL∑m=−L/2+1L/2expimφPC~mφSKmsρO/c1.The accuracy ofGcasρ,ρS in (42) depends on the number of equivalent sources L, as was shown for the case of the empty casing in Section 6.2.2, on the discretization error of Ghomρ,ρS and on M. We compared the values of Gcasρ,ρS with the exact solution as a function of N: for Muscle1 with L=M=32, the NRMSE is almost identical to that in Figure 3(a); for Leg2 with L=M=64 and for Muscle4 with L=M=128, the NRMSE is approximately 1.5 times as high as the values in Figure 3(a). When the number of equivalent sources for Muscle4 is reduced to L=M=64, the NRMSE for N=512 increases from 1.6 × 10−4 to 3.3 × 10−4, as could be expected from Table 5. Figure 8 shows the exact solution for the total field in the casing for Leg2.Figure 8
Amplitude (a, c) and phase (b, d) of the field for Leg2 in water without casing (a, b) and with casing (c, d). The field is displayed over a square subregion of width 36 cm.
(a)
(b)
(c)
(d)


#### 6.3.3. CPU Time

All summations containing complex exponentials, such as (35), (38), (43) and the right hand side in (40), are computed by means of FFTs; hence, the computational effort involved in these steps is negligible. The effort for the computation of the field on the T/R circle GcasρR,ρS is primarily determined by the solution of (40), for which we used LU factorization: we observed CPU times of 5 ms for Muscle1 with M=32, 30 ms for Leg2 with M=64, and 0.2 s for Muscle4 with M=128. The effort needed for the summations in (42) to compute the field on the mesh Gcasρ,ρS is much more important but remains lower than that needed for the CGFFT solution, as can be seen from the circles in Figure 4. It can be concluded that the embedding approach does not increase the CPU time with respect to the forward problem solution in a homogeneous environment, when only the fields on the observation contour are computed, and that it leads to a relatively small increase when the fields on the mesh are computed as well. The relative increase is less than 10% for the previously specified examples Muscle1 and Leg2 and less than 50% for Muscle 4.


## 6.1. Homogeneous Environment

### 6.1.1. Discretization Error—Choice ofN

We compared the total fieldGhomρ,ρS in the grid points, computed by solving (5) with the CGFFT method for the source position φS=0, with the exact solution as a function of the discretization cell size h. For the nonsymmetrical Leg2 object, the exact solution was chosen as the discretized solution for a much smaller cell size. We also compared the corresponding scattered field GhomscatρR,ρS in the 64 receiver points with the exact solution. Figure 3(a) shows the discretization error on the total field, with h ranging from 0.0352 cm≈λ1/224 to 2.25 cm ≈λ1/3.5. Depending on the diameter of the object, we used a mesh with side 2ℓ= 9 cm, 18 cm, or 36 cm, and we varied N between 4 and 512 (or Ngrid between 25 and 263,169). Figure 3(b) shows the resulting discretization error on the scattered field. It is clearly illustrated that the discretization errors are of Oh2 as h↓0, and it furthermore appears that for fixed h, the variation in the errors for the different cylinders considered is limited. Let us add here that for highly contrasting objects, such as air cylinders, extremely small cells may be needed in areas where the field rapidly varies—this occurs, for example, when the air interface is in the neighborhood of the source. For objects with characteristics within the ranges considered in Table 1, Figure 3(b) is helpful for choosing the cell size corresponding to a given acceptable error on the scattered field, and Figure 3(a) then gives an indication for the CGFFT stop criterion. Proceeding with the iterations when the CGFFT NRMSE (with respect to the incident field on the mesh) is smaller than a certain fraction of the total field discretization error is a waste of computational effort.Figure 3
The NRMSE for (a)Ghomρ,ρS on the mesh and for (b) GhomscatρR,ρS on ∂DO as a function of the cell size h for the objects of Table 1. In presence of the casing and for proper choices of M and L, the values in (a) are multiplied with a factor between 1 and 2 for the NRMSE of Gcasρ,ρS, and (b) remains valid for the NRMSE of GcasdifρR,ρS.
(a)
(b)


### 6.1.2. CPU Time

The CPU times for the different steps in the forward problem are illustrated in Figure4 for Leg2 as a function of N. Most of the time is taken up by the CGFFT iterations to solve (5) for the different angles of incidence; only a few percent of this time is needed to compute, respectively, the initial estimates for the total field by means of the “marching on in angle” procedure, the Green functions K in (5), and the K scattered fields GhomscatρR,ρS in the K receivers. The second most time-consuming step is the computation of the K incident fields G1ρ,ρS on the grid; the values presented here correspond to a worst-case scenario, since we did not exploit the symmetries in the grid and T/R points. Note that the Green functions can be stored for a given configuration of the mesh, transmitters/receivers, and exterior medium permittivity.Figure 4
CPU times for the different steps in the forward problem for Leg2 withK=64 sources: (1) solution of equation (5) for Ghomρ,ρS with CGFFT and “marching on in angle”; (2) computation of the incident field G1ρ,ρS; (3) computation of the scattered field GhomscatρR,ρS; (4) one CGFFT iteration; (5) the summation (42) to obtain the total field Gcasρ,ρS. The curve (6) gives the CPU time for step (1) with the incident field as initial guess.


### 6.1.3. Marching on in Angle

We tested the efficiency of “marching on in angle,” where we used the total field solutions from the three previous excitations to compute the initial estimate, by comparing the total number of CGFFT iterations to those in a conventional CGFFT forward problem solution, where the initial estimate is chosen equal to the incident field. Such a study has not yet been reported for the case of a relatively large value of the exterior medium permittivityε1r. When the contrast εrρ−ε1r/ε1r is small, only a limited number of CGFFT iterations is needed to obtain convergence.In all our tests, “marching on in angle” was faster, except for the case where the spacing between the sources was larger thanλ1. With proper choices for the number of sources K, for N and for the CGFFT stop criterium, “marching on in angle” leads to a significant reduction in the computational effort. For example, with N corresponding to a scattered field discretization error of a few percent (see Figure 3(b)), the relative reduction typically is 20% for Muscle1 N=16,K=32, 52% for Leg2 N=32,K=64—see also Figure 4 where the squares are below the dotted curve—and 32% for Muscle4 N=64,K=128, see Table 2, which gives the relative reduction in the total number of CGFFT iterations as a function of N and K. The gain in efficiency is not as spectacular as with some examples in [13], for which a much larger number of CGFFT iterations is needed to obtain convergence, which leaves more room for improving the computational efficiency.Table 2
Comparison of the total number of CGFFT iterations as a function ofN and K for the conventional (i.e., incident field initial guess) and “marching on in angle” approaches. The relative reduction in the number of iterations, the CGFFT stop criterium (NMRS CGFFT), and the resulting NRMSE for Ghom on the grid also are indicated.
N
16
32
64
128
256
512
Muscle1
K = 16
No marching
92
108
112
136
Marching
96
111
127
139
Reduction (%)
−4
—
—
—
K = 32
No marching
188
220
224
280
312
344
Marching
151
163
216
256
267
320
Reduction (%)
20
26
4
9
14
7
K = 64
No marching
380
444
448
568
632
696
Marching
176
200
325
390
460
518
Reduction (%)
54
55
27
31
27
26
NRMSE CGFFT
2 × 10−3
10−3
2 × 10−4
5 × 10−5
10−5
2 × 10−6
NRMSEGhom
5.9 × 10−3
1.7 × 10−3
4.4 × 10−4
1.1 × 10−4
2.8 × 10−5
7.1 × 10−6
Leg2
K = 32
No marching
1043
1272
1571
Marching
780
1015
1406
Reduction (%)
25
20
11
K = 64
No marching
2095
2555
3170
3489
3755
3981
Marching
878
1216
1896
2673
3239
3607
Reduction (%)
58
52
40
23
14
9
K = 128
No marching
4199
5112
6352
7533
Marching
1856
2466
3135
4596
Reduction (%)
56
52
49
39
NRMSE CGFFT
10−2
5 × 10−3
10−3
2 × 10−4
5 × 10−5
10−5
NRMSEGhom
1.1 × 10−1
3.2 × 10−2
7.9 × 10−3
1.9 × 10−3
5.1 × 10−4
1.2 × 10−4
Muscle4
K = 64
No marching
700
832
1128
1376
1648
1752
Marching
619
761
1032
1302
1563
1678
Reduction (%)
12
9
9
10
5
4
K = 128
No marching
1404
1664
2256
2760
3312
Marching
791
911
1531
2000
2586
Reduction (%)
44
45
32
28
22
K = 256
No marching
2804
3328
4512
5536
Marching
1034
1103
1679
2565
Reduction (%)
63
67
63
54
NRMSE CGFFT
2 × 10−2
10−2
2 × 10−3
5 × 10−4
10−4
5 × 10−5
NRMSEGhom
6.4 × 10−2
2.6 × 10−2
6.7 × 10−3
1.7 × 10−3
4.5 × 10−4
1.3 × 10−4It appears that the relative reduction decreases with increasing values ofN, that is, finer meshing. This can be explained as follows. The dimension of the vector space in which the solution needs to be determined increases as ON2, while the discretization error, hence the recommended CGFFT stop criterion, decreases as ON−2. The error on the “marching on in angle” initial estimate does not depend on N, while the number of iterations, which is needed to further reduce the error, increases and becomes less dependent on the choice of the initial guess.It can furthermore be seen that the relative reduction increases with increasing numbers of sourcesK. For example, when K is doubled in the aforementioned examples, the relative reduction typically is 54% for Muscle1, 52% for Leg2, and 63% for Muscle4. The sources are closer, such that the error in the “marching on in angle” initial estimate is reduced. This is also illustrated in Figure 5, which shows the number of CGFFT iterations, with and without “marching on in angle” as a function of K for Muscle1.Figure 5
The total number of CGFFT iterations as a function of the number of sourcesK for the “marching on in angle” and conventional approaches, for Muscle1 with N=16.Let us conclude by stressing the importance of choosing the largest possible value for the CGFFT stop criterion (NRMSE CGFFT), in order to get the most benefit from “marching on in angle.” When, for example, for Leg2N=32,K=64, the CGFFT stop criterion is reduced by a factor of 10, the number of CGFFT iterations is almost doubled, from 1216 to 2148, while the resulting reduction in the NRMSE on the total field, from 3.2 to 2.8%, is not significant.


### 6.1.4. Scattering Matrix—Choice ofM

For the computation of the scattering matrix, (10) is expressed in M transmitter/receiver positions on ∂DO, and the summations are truncated to M terms,
(31)GhomscatρR,ρS=12π∑m=−M/2+1M/2∑m′=−M/2+1M/2Sm,m′KmsρOc1Km′sρOc1expi2πmr−m′tM,with r,t=1,…,M, and where M should be large enough for the aliasing effects to be negligible or Sm+lM,m′+l′MKm+lMKm′+l′M≪Sm,m′KmKm′ for ±l,l′=1,…,∞. The elements
(32)Qm,m′=M2πSm,m′KmsρOc1Km′sρOc1are computed from (31) by means of a 2D FFT. The computational effort of this step is negligible: 10 ms for the example of Figure 4. For a centered circular homogeneous cylinder, the scattering matrix is diagonal with Qm,m given by
(33)−M2πc1Insρ2/c1In−1sρ2/c2−c2Insρ2/c2In−1sρ2/c1c1Knsρ2/c1In−1sρ2/c2−c2Insρ2/c2Kn−1sρ2/c1Km2sρOc1,with ρ2 and c2=c0/ϵ2r the radius of and the wave speed in the cylinder, respectively. For high orders m, the elements (33) decrease with increasing order at a rate which primarily depends on ρ2, for a given ρO. This is illustrated in Table 3, which shows that the truncation errors in (31) are negligible when M=32>2ODP for the smallest cylinders, which yields double precision, and M=64>2OSP for the largest cylinder, which is largely sufficient to obtain single precision. The accuracy of the scattering matrix then is determined by the discretization error on Ghomscat.Table 3
The ordersm=OSP and m=ODP for which the elements Qmm of the scattering matrix have dropped to fractions 10−7 (single precision) and 10−14 (double precision), respectively, of their maximum absolute value, for the homogeneous cylinders of Table 1.
Cylinder
Muscle1
Air1
Air2
Muscle4
OSP
9
9
14
24
ODP
14
14
20
36


## 6.1.1. Discretization Error—Choice ofN

We compared the total fieldGhomρ,ρS in the grid points, computed by solving (5) with the CGFFT method for the source position φS=0, with the exact solution as a function of the discretization cell size h. For the nonsymmetrical Leg2 object, the exact solution was chosen as the discretized solution for a much smaller cell size. We also compared the corresponding scattered field GhomscatρR,ρS in the 64 receiver points with the exact solution. Figure 3(a) shows the discretization error on the total field, with h ranging from 0.0352 cm≈λ1/224 to 2.25 cm ≈λ1/3.5. Depending on the diameter of the object, we used a mesh with side 2ℓ= 9 cm, 18 cm, or 36 cm, and we varied N between 4 and 512 (or Ngrid between 25 and 263,169). Figure 3(b) shows the resulting discretization error on the scattered field. It is clearly illustrated that the discretization errors are of Oh2 as h↓0, and it furthermore appears that for fixed h, the variation in the errors for the different cylinders considered is limited. Let us add here that for highly contrasting objects, such as air cylinders, extremely small cells may be needed in areas where the field rapidly varies—this occurs, for example, when the air interface is in the neighborhood of the source. For objects with characteristics within the ranges considered in Table 1, Figure 3(b) is helpful for choosing the cell size corresponding to a given acceptable error on the scattered field, and Figure 3(a) then gives an indication for the CGFFT stop criterion. Proceeding with the iterations when the CGFFT NRMSE (with respect to the incident field on the mesh) is smaller than a certain fraction of the total field discretization error is a waste of computational effort.Figure 3
The NRMSE for (a)Ghomρ,ρS on the mesh and for (b) GhomscatρR,ρS on ∂DO as a function of the cell size h for the objects of Table 1. In presence of the casing and for proper choices of M and L, the values in (a) are multiplied with a factor between 1 and 2 for the NRMSE of Gcasρ,ρS, and (b) remains valid for the NRMSE of GcasdifρR,ρS.
(a)
(b)


## 6.1.2. CPU Time

The CPU times for the different steps in the forward problem are illustrated in Figure4 for Leg2 as a function of N. Most of the time is taken up by the CGFFT iterations to solve (5) for the different angles of incidence; only a few percent of this time is needed to compute, respectively, the initial estimates for the total field by means of the “marching on in angle” procedure, the Green functions K in (5), and the K scattered fields GhomscatρR,ρS in the K receivers. The second most time-consuming step is the computation of the K incident fields G1ρ,ρS on the grid; the values presented here correspond to a worst-case scenario, since we did not exploit the symmetries in the grid and T/R points. Note that the Green functions can be stored for a given configuration of the mesh, transmitters/receivers, and exterior medium permittivity.Figure 4
CPU times for the different steps in the forward problem for Leg2 withK=64 sources: (1) solution of equation (5) for Ghomρ,ρS with CGFFT and “marching on in angle”; (2) computation of the incident field G1ρ,ρS; (3) computation of the scattered field GhomscatρR,ρS; (4) one CGFFT iteration; (5) the summation (42) to obtain the total field Gcasρ,ρS. The curve (6) gives the CPU time for step (1) with the incident field as initial guess.


## 6.1.3. Marching on in Angle

We tested the efficiency of “marching on in angle,” where we used the total field solutions from the three previous excitations to compute the initial estimate, by comparing the total number of CGFFT iterations to those in a conventional CGFFT forward problem solution, where the initial estimate is chosen equal to the incident field. Such a study has not yet been reported for the case of a relatively large value of the exterior medium permittivityε1r. When the contrast εrρ−ε1r/ε1r is small, only a limited number of CGFFT iterations is needed to obtain convergence.In all our tests, “marching on in angle” was faster, except for the case where the spacing between the sources was larger thanλ1. With proper choices for the number of sources K, for N and for the CGFFT stop criterium, “marching on in angle” leads to a significant reduction in the computational effort. For example, with N corresponding to a scattered field discretization error of a few percent (see Figure 3(b)), the relative reduction typically is 20% for Muscle1 N=16,K=32, 52% for Leg2 N=32,K=64—see also Figure 4 where the squares are below the dotted curve—and 32% for Muscle4 N=64,K=128, see Table 2, which gives the relative reduction in the total number of CGFFT iterations as a function of N and K. The gain in efficiency is not as spectacular as with some examples in [13], for which a much larger number of CGFFT iterations is needed to obtain convergence, which leaves more room for improving the computational efficiency.Table 2
Comparison of the total number of CGFFT iterations as a function ofN and K for the conventional (i.e., incident field initial guess) and “marching on in angle” approaches. The relative reduction in the number of iterations, the CGFFT stop criterium (NMRS CGFFT), and the resulting NRMSE for Ghom on the grid also are indicated.
N
16
32
64
128
256
512
Muscle1
K = 16
No marching
92
108
112
136
Marching
96
111
127
139
Reduction (%)
−4
—
—
—
K = 32
No marching
188
220
224
280
312
344
Marching
151
163
216
256
267
320
Reduction (%)
20
26
4
9
14
7
K = 64
No marching
380
444
448
568
632
696
Marching
176
200
325
390
460
518
Reduction (%)
54
55
27
31
27
26
NRMSE CGFFT
2 × 10−3
10−3
2 × 10−4
5 × 10−5
10−5
2 × 10−6
NRMSEGhom
5.9 × 10−3
1.7 × 10−3
4.4 × 10−4
1.1 × 10−4
2.8 × 10−5
7.1 × 10−6
Leg2
K = 32
No marching
1043
1272
1571
Marching
780
1015
1406
Reduction (%)
25
20
11
K = 64
No marching
2095
2555
3170
3489
3755
3981
Marching
878
1216
1896
2673
3239
3607
Reduction (%)
58
52
40
23
14
9
K = 128
No marching
4199
5112
6352
7533
Marching
1856
2466
3135
4596
Reduction (%)
56
52
49
39
NRMSE CGFFT
10−2
5 × 10−3
10−3
2 × 10−4
5 × 10−5
10−5
NRMSEGhom
1.1 × 10−1
3.2 × 10−2
7.9 × 10−3
1.9 × 10−3
5.1 × 10−4
1.2 × 10−4
Muscle4
K = 64
No marching
700
832
1128
1376
1648
1752
Marching
619
761
1032
1302
1563
1678
Reduction (%)
12
9
9
10
5
4
K = 128
No marching
1404
1664
2256
2760
3312
Marching
791
911
1531
2000
2586
Reduction (%)
44
45
32
28
22
K = 256
No marching
2804
3328
4512
5536
Marching
1034
1103
1679
2565
Reduction (%)
63
67
63
54
NRMSE CGFFT
2 × 10−2
10−2
2 × 10−3
5 × 10−4
10−4
5 × 10−5
NRMSEGhom
6.4 × 10−2
2.6 × 10−2
6.7 × 10−3
1.7 × 10−3
4.5 × 10−4
1.3 × 10−4It appears that the relative reduction decreases with increasing values ofN, that is, finer meshing. This can be explained as follows. The dimension of the vector space in which the solution needs to be determined increases as ON2, while the discretization error, hence the recommended CGFFT stop criterion, decreases as ON−2. The error on the “marching on in angle” initial estimate does not depend on N, while the number of iterations, which is needed to further reduce the error, increases and becomes less dependent on the choice of the initial guess.It can furthermore be seen that the relative reduction increases with increasing numbers of sourcesK. For example, when K is doubled in the aforementioned examples, the relative reduction typically is 54% for Muscle1, 52% for Leg2, and 63% for Muscle4. The sources are closer, such that the error in the “marching on in angle” initial estimate is reduced. This is also illustrated in Figure 5, which shows the number of CGFFT iterations, with and without “marching on in angle” as a function of K for Muscle1.Figure 5
The total number of CGFFT iterations as a function of the number of sourcesK for the “marching on in angle” and conventional approaches, for Muscle1 with N=16.Let us conclude by stressing the importance of choosing the largest possible value for the CGFFT stop criterion (NRMSE CGFFT), in order to get the most benefit from “marching on in angle.” When, for example, for Leg2N=32,K=64, the CGFFT stop criterion is reduced by a factor of 10, the number of CGFFT iterations is almost doubled, from 1216 to 2148, while the resulting reduction in the NRMSE on the total field, from 3.2 to 2.8%, is not significant.


## 6.1.4. Scattering Matrix—Choice ofM

For the computation of the scattering matrix, (10) is expressed in M transmitter/receiver positions on ∂DO, and the summations are truncated to M terms,
(31)GhomscatρR,ρS=12π∑m=−M/2+1M/2∑m′=−M/2+1M/2Sm,m′KmsρOc1Km′sρOc1expi2πmr−m′tM,with r,t=1,…,M, and where M should be large enough for the aliasing effects to be negligible or Sm+lM,m′+l′MKm+lMKm′+l′M≪Sm,m′KmKm′ for ±l,l′=1,…,∞. The elements
(32)Qm,m′=M2πSm,m′KmsρOc1Km′sρOc1are computed from (31) by means of a 2D FFT. The computational effort of this step is negligible: 10 ms for the example of Figure 4. For a centered circular homogeneous cylinder, the scattering matrix is diagonal with Qm,m given by
(33)−M2πc1Insρ2/c1In−1sρ2/c2−c2Insρ2/c2In−1sρ2/c1c1Knsρ2/c1In−1sρ2/c2−c2Insρ2/c2Kn−1sρ2/c1Km2sρOc1,with ρ2 and c2=c0/ϵ2r the radius of and the wave speed in the cylinder, respectively. For high orders m, the elements (33) decrease with increasing order at a rate which primarily depends on ρ2, for a given ρO. This is illustrated in Table 3, which shows that the truncation errors in (31) are negligible when M=32>2ODP for the smallest cylinders, which yields double precision, and M=64>2OSP for the largest cylinder, which is largely sufficient to obtain single precision. The accuracy of the scattering matrix then is determined by the discretization error on Ghomscat.Table 3
The ordersm=OSP and m=ODP for which the elements Qmm of the scattering matrix have dropped to fractions 10−7 (single precision) and 10−14 (double precision), respectively, of their maximum absolute value, for the homogeneous cylinders of Table 1.
Cylinder
Muscle1
Air1
Air2
Muscle4
OSP
9
9
14
24
ODP
14
14
20
36


## 6.2. Direct Reflection by the Casing

The direct reflection from the casing, which is treated as a separate constituent in the expressions for the field on∂DO and for the equivalent surface current, does not depend on the object. Hence, these contributions can be computed once beforehand for a given water-filled scanner geometry. Let us therefore consider the case of a line source on ∂DO in the casing without object. The spectral representation of the field Gcasemptyρ,ρS in the complete region 0<ρ<b is then given by the first two sums in (26).


### 6.2.1. Exact Solution

The convergence of the first sum in (26) is extremely slow for observation points on ∂DO. This is visible in Figure 6, in Table 4, and from applying the asymptotic approximation ([34], 9.3.1) for high orders m≫∣s/c1ρ>∣,
(34)12πImsρc1KmsρOc1≈14πmρ<ρ>m.Figure 6
Absolute values of the terms in (26) for ρ=ρO: (1) 1/2π∣ImKm∣, (2) 1/2π∣RmImIm∣, ∣C~mIm∣ computed with equation (40) for (3) Muscle1, (4) Leg2, and (5) Muscle4. The terms ∣DmKm∣ are indicated with +.Table 4
The ordersm=OSP and m=ODP for which the terms in the series in (8) and (35) have dropped to fractions 10-7 (single precision) and 10-14 (double precision), respectively, of their maximum absolute values, for different ρ.
ρ (cm)
27.6
25
15
5
(8): ∣Imsρ/c1KmsρO/c1∣
ρ/ρO
1
0.91
0.54
0.18
OSP
≈107
130
30
15
ODP
≈1014
283
52
22
(35): ∣RmImsρ/c1ImsρO/c1∣
ρ/ρ1
0.91
0.82
0.49
0.16
OSP
97
60
27
15
ODP
211
124
44
22Consequently, it is more efficient to computeG1ρ,ρS with (3) or
(35)Gcasemptyρ,ρS=12πK0sc1ρ−ρS+12π∑m=−OOexpimφ−φSRmImsρOc1Imsρc1,where the reflected-field constituent is truncated at order O. The convergence of the spectral representation of the direct casing reflection in (35) is better than that of (8) in the region 0<ρ<ρO; see Figure 6 and Table 4. In this case, the asymptotic approximation of the term m, for m≫∣s/c1b∣, yields
(36)12πRmImsρOc1Imsρc1→−14πmρρOb2m,which is identical to that of a line source located on a contour with radius ρ1=b2/ρO in a homogeneous background—b thus is the geometric mean of ρO and ρ1. We observed that the direct casing reflection reaches stable and accurate DP values in the 64 receivers on ∂DO when choosing O=ODP=211.


### 6.2.2. Expression from Equivalent Line Sources

Alternatively, in the regionρ<ρO, the field can be regarded as if it originated from an equivalent surface current (21), which we replace with a discrete set of L equivalent line sources with spacing ΔφP=2π/L on ∂DO in a homogeneous background, denoting by WemptyφP,φS the complex amplitude of the equivalent source at φP due to an excitation at φS. Imposing the identity of the reflected field constituent in (35) and the field generated by these sources, for ρ<ρO,
(37)12π∑m=−L/2+1L/2expimφ−φSRmImsρOc1Imsρc1=12π∑φPWemptyφP,φS∑m=−L/2+1L/2expimφ−φp⋅Imsρc1KmsρOc1,where we used the spectral representation (8) for G1ρ,ρP, we obtain the following expression for the complex amplitude
(38)WemptyφP,φS=1L∑m=−L/2+1L/2expimφP−φSRmImsρO/c1KmsρO/c1.According to (23) and since in this section Ghomρ,ρS=G1ρ,ρS, the field in the region ρ<ρO then is expressed as the linear combination
(39)Gcasemptyρ,ρS=G1ρ,ρS+∑φPWemptyφP,φSG1ρ,ρP.The truncation toL terms in the RHS of (37) may lead to significant errors on the field (39) in grid points that are near to ∂DO. In Table 5, we compared the field on the mesh (39), generated by different numbers L of equivalent line sources, with the exact solution (35) for meshes with sides 9 cm, 18 cm, and 36 cm. It follows that for the mesh with side 9 cm, L=32 already yields a very high precision; for the mesh with side 36 cm, L=64 yields a moderate precision, and L≥128 is needed for high precision computations. Figure 7 shows an image of the amplitude and phase of the field (35) in the water-filled casing.Table 5
The NRMSE of the field (39) in the empty casing as a function of the number of equivalent sources L for different sizes 2ℓ of the mesh.
2ℓ (cm)
L=16
L=32
L=64
L=128
L=256
9
5 ×10−4
2 ×10−9
2 ×10−14
2 ×10−14
18
6 ×10−2
7 ×10−5
4 ×10−14
4 ×10−14
36
10−1
2 ×10−4
3 ×10−7
5 ×10−11Figure 7
Amplitude (a) and phase (b) of the field (35) in the empty water-filled casing. An upper limit of 0.2 was chosen in the image of the amplitude.
(a)
(b)


## 6.2.1. Exact Solution

The convergence of the first sum in (26) is extremely slow for observation points on ∂DO. This is visible in Figure 6, in Table 4, and from applying the asymptotic approximation ([34], 9.3.1) for high orders m≫∣s/c1ρ>∣,
(34)12πImsρc1KmsρOc1≈14πmρ<ρ>m.Figure 6
Absolute values of the terms in (26) for ρ=ρO: (1) 1/2π∣ImKm∣, (2) 1/2π∣RmImIm∣, ∣C~mIm∣ computed with equation (40) for (3) Muscle1, (4) Leg2, and (5) Muscle4. The terms ∣DmKm∣ are indicated with +.Table 4
The ordersm=OSP and m=ODP for which the terms in the series in (8) and (35) have dropped to fractions 10-7 (single precision) and 10-14 (double precision), respectively, of their maximum absolute values, for different ρ.
ρ (cm)
27.6
25
15
5
(8): ∣Imsρ/c1KmsρO/c1∣
ρ/ρO
1
0.91
0.54
0.18
OSP
≈107
130
30
15
ODP
≈1014
283
52
22
(35): ∣RmImsρ/c1ImsρO/c1∣
ρ/ρ1
0.91
0.82
0.49
0.16
OSP
97
60
27
15
ODP
211
124
44
22Consequently, it is more efficient to computeG1ρ,ρS with (3) or
(35)Gcasemptyρ,ρS=12πK0sc1ρ−ρS+12π∑m=−OOexpimφ−φSRmImsρOc1Imsρc1,where the reflected-field constituent is truncated at order O. The convergence of the spectral representation of the direct casing reflection in (35) is better than that of (8) in the region 0<ρ<ρO; see Figure 6 and Table 4. In this case, the asymptotic approximation of the term m, for m≫∣s/c1b∣, yields
(36)12πRmImsρOc1Imsρc1→−14πmρρOb2m,which is identical to that of a line source located on a contour with radius ρ1=b2/ρO in a homogeneous background—b thus is the geometric mean of ρO and ρ1. We observed that the direct casing reflection reaches stable and accurate DP values in the 64 receivers on ∂DO when choosing O=ODP=211.


## 6.2.2. Expression from Equivalent Line Sources

Alternatively, in the regionρ<ρO, the field can be regarded as if it originated from an equivalent surface current (21), which we replace with a discrete set of L equivalent line sources with spacing ΔφP=2π/L on ∂DO in a homogeneous background, denoting by WemptyφP,φS the complex amplitude of the equivalent source at φP due to an excitation at φS. Imposing the identity of the reflected field constituent in (35) and the field generated by these sources, for ρ<ρO,
(37)12π∑m=−L/2+1L/2expimφ−φSRmImsρOc1Imsρc1=12π∑φPWemptyφP,φS∑m=−L/2+1L/2expimφ−φp⋅Imsρc1KmsρOc1,where we used the spectral representation (8) for G1ρ,ρP, we obtain the following expression for the complex amplitude
(38)WemptyφP,φS=1L∑m=−L/2+1L/2expimφP−φSRmImsρO/c1KmsρO/c1.According to (23) and since in this section Ghomρ,ρS=G1ρ,ρS, the field in the region ρ<ρO then is expressed as the linear combination
(39)Gcasemptyρ,ρS=G1ρ,ρS+∑φPWemptyφP,φSG1ρ,ρP.The truncation toL terms in the RHS of (37) may lead to significant errors on the field (39) in grid points that are near to ∂DO. In Table 5, we compared the field on the mesh (39), generated by different numbers L of equivalent line sources, with the exact solution (35) for meshes with sides 9 cm, 18 cm, and 36 cm. It follows that for the mesh with side 9 cm, L=32 already yields a very high precision; for the mesh with side 36 cm, L=64 yields a moderate precision, and L≥128 is needed for high precision computations. Figure 7 shows an image of the amplitude and phase of the field (35) in the water-filled casing.Table 5
The NRMSE of the field (39) in the empty casing as a function of the number of equivalent sources L for different sizes 2ℓ of the mesh.
2ℓ (cm)
L=16
L=32
L=64
L=128
L=256
9
5 ×10−4
2 ×10−9
2 ×10−14
2 ×10−14
18
6 ×10−2
7 ×10−5
4 ×10−14
4 ×10−14
36
10−1
2 ×10−4
3 ×10−7
5 ×10−11Figure 7
Amplitude (a) and phase (b) of the field (35) in the empty water-filled casing. An upper limit of 0.2 was chosen in the image of the amplitude.
(a)
(b)


## 6.3. Complete Configuration

Finally, we look into the computation of the fields for the complete configuration of an object in the casing. It is shown that the embedding approach maintains the accuracy of the forward problem solution in homogeneous space, ifM (order of the scattering matrix) and L (number of equivalent line sources) are properly chosen. The additional computational effort for the embedding as such is also examined.


### 6.3.1. Scattered Field on the Receivers

Instead of solving (16) for Cm, we take into account the separation of the direct casing reflection (26) and solve the resulting set of equations, truncated to M terms, for C~mIm
(40)C~mφSImsρOc1−Um∑n=−M/2+1M/22πVnQmnMC~nφSInsρOc1=Um∑n=−M/2+1M/2QmnM1+Unexp−inφS,where Um=RmIm/Km and Vm=ImKm with m=−M/2+1…M/2. The elements DmKm then are computed with
(41)DmφSKmsρOc1=C~mφSImsρO/c1Um.The absolute values ofC~mIm and DmKm rapidly decrease as a function of m, as is shown in Figure 6. As a consequence, the convergence of the spectral representation (26) of the field GcasdifρR,ρS is much better than that of GcasemptyρR,ρS. Provided that M is chosen according to Table 3, the error on GcasdifρR,ρS only depends on the discretization error in the scattered field GhomscatρR,ρS. We compared the values of GcasdifρR,ρS with the exact solution as a function of N for Muscle1 with M=32, Leg2 with M=64, and Muscle4 with M=128, and we observed that the resulting NRMSE is almost identical to that for the homogeneous solutions; see Figure 3(b).


### 6.3.2. Total Field on the Grid

The total fieldGcasρ,ρS on the mesh is computed as a linear combination of homogeneous solutions for L equivalent line sources
(42)Gcasρ,ρS=Ghomρ,ρS+∑φPWemptyφP,φS+WdifφP,φSGhomρ,ρP,where WemptyφP,φS is given by (38) and where
(43)WdifφP,φS=2πL∑m=−L/2+1L/2expimφPC~mφSKmsρO/c1.The accuracy ofGcasρ,ρS in (42) depends on the number of equivalent sources L, as was shown for the case of the empty casing in Section 6.2.2, on the discretization error of Ghomρ,ρS and on M. We compared the values of Gcasρ,ρS with the exact solution as a function of N: for Muscle1 with L=M=32, the NRMSE is almost identical to that in Figure 3(a); for Leg2 with L=M=64 and for Muscle4 with L=M=128, the NRMSE is approximately 1.5 times as high as the values in Figure 3(a). When the number of equivalent sources for Muscle4 is reduced to L=M=64, the NRMSE for N=512 increases from 1.6 × 10−4 to 3.3 × 10−4, as could be expected from Table 5. Figure 8 shows the exact solution for the total field in the casing for Leg2.Figure 8
Amplitude (a, c) and phase (b, d) of the field for Leg2 in water without casing (a, b) and with casing (c, d). The field is displayed over a square subregion of width 36 cm.
(a)
(b)
(c)
(d)


### 6.3.3. CPU Time

All summations containing complex exponentials, such as (35), (38), (43) and the right hand side in (40), are computed by means of FFTs; hence, the computational effort involved in these steps is negligible. The effort for the computation of the field on the T/R circle GcasρR,ρS is primarily determined by the solution of (40), for which we used LU factorization: we observed CPU times of 5 ms for Muscle1 with M=32, 30 ms for Leg2 with M=64, and 0.2 s for Muscle4 with M=128. The effort needed for the summations in (42) to compute the field on the mesh Gcasρ,ρS is much more important but remains lower than that needed for the CGFFT solution, as can be seen from the circles in Figure 4. It can be concluded that the embedding approach does not increase the CPU time with respect to the forward problem solution in a homogeneous environment, when only the fields on the observation contour are computed, and that it leads to a relatively small increase when the fields on the mesh are computed as well. The relative increase is less than 10% for the previously specified examples Muscle1 and Leg2 and less than 50% for Muscle 4.


## 6.3.1. Scattered Field on the Receivers

Instead of solving (16) for Cm, we take into account the separation of the direct casing reflection (26) and solve the resulting set of equations, truncated to M terms, for C~mIm
(40)C~mφSImsρOc1−Um∑n=−M/2+1M/22πVnQmnMC~nφSInsρOc1=Um∑n=−M/2+1M/2QmnM1+Unexp−inφS,where Um=RmIm/Km and Vm=ImKm with m=−M/2+1…M/2. The elements DmKm then are computed with
(41)DmφSKmsρOc1=C~mφSImsρO/c1Um.The absolute values ofC~mIm and DmKm rapidly decrease as a function of m, as is shown in Figure 6. As a consequence, the convergence of the spectral representation (26) of the field GcasdifρR,ρS is much better than that of GcasemptyρR,ρS. Provided that M is chosen according to Table 3, the error on GcasdifρR,ρS only depends on the discretization error in the scattered field GhomscatρR,ρS. We compared the values of GcasdifρR,ρS with the exact solution as a function of N for Muscle1 with M=32, Leg2 with M=64, and Muscle4 with M=128, and we observed that the resulting NRMSE is almost identical to that for the homogeneous solutions; see Figure 3(b).


## 6.3.2. Total Field on the Grid

The total fieldGcasρ,ρS on the mesh is computed as a linear combination of homogeneous solutions for L equivalent line sources
(42)Gcasρ,ρS=Ghomρ,ρS+∑φPWemptyφP,φS+WdifφP,φSGhomρ,ρP,where WemptyφP,φS is given by (38) and where
(43)WdifφP,φS=2πL∑m=−L/2+1L/2expimφPC~mφSKmsρO/c1.The accuracy ofGcasρ,ρS in (42) depends on the number of equivalent sources L, as was shown for the case of the empty casing in Section 6.2.2, on the discretization error of Ghomρ,ρS and on M. We compared the values of Gcasρ,ρS with the exact solution as a function of N: for Muscle1 with L=M=32, the NRMSE is almost identical to that in Figure 3(a); for Leg2 with L=M=64 and for Muscle4 with L=M=128, the NRMSE is approximately 1.5 times as high as the values in Figure 3(a). When the number of equivalent sources for Muscle4 is reduced to L=M=64, the NRMSE for N=512 increases from 1.6 × 10−4 to 3.3 × 10−4, as could be expected from Table 5. Figure 8 shows the exact solution for the total field in the casing for Leg2.Figure 8
Amplitude (a, c) and phase (b, d) of the field for Leg2 in water without casing (a, b) and with casing (c, d). The field is displayed over a square subregion of width 36 cm.
(a)
(b)
(c)
(d)


## 6.3.3. CPU Time

All summations containing complex exponentials, such as (35), (38), (43) and the right hand side in (40), are computed by means of FFTs; hence, the computational effort involved in these steps is negligible. The effort for the computation of the field on the T/R circle GcasρR,ρS is primarily determined by the solution of (40), for which we used LU factorization: we observed CPU times of 5 ms for Muscle1 with M=32, 30 ms for Leg2 with M=64, and 0.2 s for Muscle4 with M=128. The effort needed for the summations in (42) to compute the field on the mesh Gcasρ,ρS is much more important but remains lower than that needed for the CGFFT solution, as can be seen from the circles in Figure 4. It can be concluded that the embedding approach does not increase the CPU time with respect to the forward problem solution in a homogeneous environment, when only the fields on the observation contour are computed, and that it leads to a relatively small increase when the fields on the mesh are computed as well. The relative increase is less than 10% for the previously specified examples Muscle1 and Leg2 and less than 50% for Muscle 4.


## 7. Computational Procedure

We conclude this paper with a summary of the computational details. In analyzing the computational complexity, it should be kept in mind that the approach described in this paper was devised for use as forward scheme in inverse profiling, where an unknown configuration is reconstructed by matching the corresponding scattered field to a known measured field by linear or nonlinear optimization. For each new estimate of the configuration, the field caused byK sources on the observation contour must be determined. When the optimization converges, the successive estimates gradually approach the desired optimum.


### 7.1. Homogeneous Environment

To demonstrate the efficiency of the scheme, we compare it with a straightforward implementation of the method of moments. For the object in a homogeneous environment, the first advantage is the second-order accuracy of the space discretization. In Figure3, the conclusion from [12] that the error in the computed fields is of Oh2 for decreasing h or equivalently of ON−2 for increasing N, was confirmed for our test objects. Second, as shown in (5), the convolution structure of the continuous equation (2) was preserved. A straightforward evaluation of a matrix-vector product requires an effort of ON4 in each CG iteration step for N2 unknown field values. Replacing these multiplications by two-dimensional FFT operations reduces the computational complexity to ON2lnNper step. Third, marching on in angle reduces the number of iterations. From Table 2, an acceleration by about 50% is observed. In all cases, the computational procedure is considerably more efficient than a straightforward matrix inversion, which requires a computational effort of ON6, followed by K matrix-vector computations at an effort of OKN4. In fact, the motivation for treating only the profile parameters as independent variables during the optimization in [5, 6, 18] was the efficiency of this forward scheme.


### 7.2. Complete Configuration

For the object inside the scanner, the conventional approach requires evaluating the field due to line sources in an empty scanner forρS at N2 mesh nodes, each for observation points ρ at N2 mesh nodes. For each pair of nodes, the modified Bessel functions have to be computed for M+1 orders, which leads to a matrix fill time of OMN4. In the algorithm described in this paper, the M+1 Bessel functions are computed once for the argument sρ0/c1, and we invert the truncated version (40) of (16). Moreover, for each line source, we synthesize the actual field (42) at N2 points from K fields in a homogeneous environment, which amounts to a total effort of OK2N2. As mentioned towards the end of Section 6.3.3, only the last step leads to a relatively small increase in computation time. The proposed embedding approach thus is significantly more efficient than computing the fields in the complete configuration by means of Green’s functions of the empty casing [19, 20].


### 7.3. Guidelines

Last but not least, we enumerate the various steps of our algorithm, giving some guidelines for efficient application. The goal is to compute with accuracyA the difference field on the T/R circle and the total field on the mesh for an object with a maximum size of 2ℓ. Based on the results given in this paper, we recommend the following procedures:
(1)
Choose the smallest possible value forN, which yields an accuracy A, with the aid of Figure 3(b).
(2)
Choose the smallest possible value forL, which yields a (much) better accuracy than A on the mesh, with the aid of Table 5.
(3)
Choose the smallest possible value forM, which yields a (much) better accuracy than A for the given object size, with the aid of Table 3.
(4)
In general,L>M; hence, choose the number of forward problem solutions K=L. For convenience, we have chosen M=K in all our examples.
(5)
ComputeRm; Wempty; and, if also the total field on the T/R circle is needed, GcasemptyρR,ρS.
(6)
Compute Green’s functions in a homogeneous environment.
(7)
Choose the CGFFT stop criterion, with the aid of Figure3(a), and solve the forward problem in a homogeneous environment with “marching on in angle.”
(8)
Compute the elementsQm,m′ of the scattering matrix and the elements C~mIm and DmKm.
(9)
Compute the field on theT/R circle.
(10)
ComputeWdif and the field on the mesh.


## 7.1. Homogeneous Environment

To demonstrate the efficiency of the scheme, we compare it with a straightforward implementation of the method of moments. For the object in a homogeneous environment, the first advantage is the second-order accuracy of the space discretization. In Figure3, the conclusion from [12] that the error in the computed fields is of Oh2 for decreasing h or equivalently of ON−2 for increasing N, was confirmed for our test objects. Second, as shown in (5), the convolution structure of the continuous equation (2) was preserved. A straightforward evaluation of a matrix-vector product requires an effort of ON4 in each CG iteration step for N2 unknown field values. Replacing these multiplications by two-dimensional FFT operations reduces the computational complexity to ON2lnNper step. Third, marching on in angle reduces the number of iterations. From Table 2, an acceleration by about 50% is observed. In all cases, the computational procedure is considerably more efficient than a straightforward matrix inversion, which requires a computational effort of ON6, followed by K matrix-vector computations at an effort of OKN4. In fact, the motivation for treating only the profile parameters as independent variables during the optimization in [5, 6, 18] was the efficiency of this forward scheme.


## 7.2. Complete Configuration

For the object inside the scanner, the conventional approach requires evaluating the field due to line sources in an empty scanner forρS at N2 mesh nodes, each for observation points ρ at N2 mesh nodes. For each pair of nodes, the modified Bessel functions have to be computed for M+1 orders, which leads to a matrix fill time of OMN4. In the algorithm described in this paper, the M+1 Bessel functions are computed once for the argument sρ0/c1, and we invert the truncated version (40) of (16). Moreover, for each line source, we synthesize the actual field (42) at N2 points from K fields in a homogeneous environment, which amounts to a total effort of OK2N2. As mentioned towards the end of Section 6.3.3, only the last step leads to a relatively small increase in computation time. The proposed embedding approach thus is significantly more efficient than computing the fields in the complete configuration by means of Green’s functions of the empty casing [19, 20].


## 7.3. Guidelines

Last but not least, we enumerate the various steps of our algorithm, giving some guidelines for efficient application. The goal is to compute with accuracyA the difference field on the T/R circle and the total field on the mesh for an object with a maximum size of 2ℓ. Based on the results given in this paper, we recommend the following procedures:
(1)
Choose the smallest possible value forN, which yields an accuracy A, with the aid of Figure 3(b).
(2)
Choose the smallest possible value forL, which yields a (much) better accuracy than A on the mesh, with the aid of Table 5.
(3)
Choose the smallest possible value forM, which yields a (much) better accuracy than A for the given object size, with the aid of Table 3.
(4)
In general,L>M; hence, choose the number of forward problem solutions K=L. For convenience, we have chosen M=K in all our examples.
(5)
ComputeRm; Wempty; and, if also the total field on the T/R circle is needed, GcasemptyρR,ρS.
(6)
Compute Green’s functions in a homogeneous environment.
(7)
Choose the CGFFT stop criterion, with the aid of Figure3(a), and solve the forward problem in a homogeneous environment with “marching on in angle.”
(8)
Compute the elementsQm,m′ of the scattering matrix and the elements C~mIm and DmKm.
(9)
Compute the field on theT/R circle.
(10)
ComputeWdif and the field on the mesh.


## 8. Conclusion

In this paper, we have described a procedure to decompose the computation of electromagnetic fields in a relatively complicated configuration. The procedure allows (re)computing the field in part of the configuration, while the remaining part of the configuration and its electromagnetic response is left unchanged. By considering scattering in a homogeneous environment, the efficiency of the CGFFT procedure is exploited. Generating the initial estimate by “marching on in angle” accelerates the convergence of this procedure significantly.The procedure has been applied to a standard scanner configuration for 2D inverse profiling. In the model, we have neglected the influence of the finite length of the antennas, the mutual coupling, and the variation in the properties of the individual antennas. Previous expertise [35] has shown that such effects can be handled in the calibration of the results, which is needed anyway. Reconstruction results for the idealized configuration have already been described in [18].In the present paper, we have addressed the efficiency and accuracy of the forward algorithm and described the influence of the different tuning parameters in the algorithm. Results have been presented and discussed for canonical objects with representative values for the permittivity and the object dimension. A systematic procedure has been proposed for choosing computational parameters such as truncation limits and tolerances.



---

*Source: 3020417-2018-05-10.xml*