An Approximation Method for Electrical Impedance Tomography by Paulo J. S. Pereira B.A.Sc., The University of British Columbia, 2008 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Applied Science in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University Of British Columbia (Vancouver, Canada) August, 2008 c© Paulo J. S. Pereira 2008 Abstract Electrical impedance tomography is an imaging method with applications to geophysics and medical imaging. A new approximation is presented based on Nachman’s 2-dimensional construction for closed domains. It improves upon existing approximations by extending the range of application from resolving 2 times the surface conductivity to imaging perfect conductors and insulators. With perfect knowledge of boundary data, this approximation exactly resolves a single conductive disc embedded in a homogenous domain. The problem, how- ever, is ill-posed, and imaging performance degrades quickly as the distance from the boundary increases. The key to the approximation lies in (a) approximating Fadeev’s Green func- tion (b) pre-processing measured voltages based on a boundary-integral equation (c) solving a linearized inverse problem (d) solving a d-bar equation, and (e) scaling the resulting image based on analytical results for a disc. In the devel- opment of the approximation, a new formula for Fadeev’s Green’s function is presented in terms of the Exponential Integral function. Also, new comparisons are made between reconstructions with and without solving the d-bar equation, showing that the added computational expense of solving the d-bar equation is not justified for radial problems. There is no discernible improvement in image quality. As a result, the approximation converts the inverse conductivity prob- lem into a novel one-step linear problem with pre-conditioning of boundary data and scaling of the resulting image. Several extensions to this work are possible. The approximation is imple- mented for a circular domain with unit conductivity near the boundary, and extensions to other domains, bounded and unbounded should be possible, with non-constant conductivity near the boundary requiring further approximation. Ultimately, further research is required to ascertain whether it is possible to extend these techniques to imaging problems in three dimensions. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v I Thesis 1 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1 Overview of Nachman’s Construction . . . . . . . . . . . . . . . 3 1.2 Organization of Work . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Survey of Electrical Impedance Tomography . . . . . . . . . . 5 2.1 Apparent Resistivity . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Variational Methods . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Resolution and Optimal Current Patterns . . . . . . . . . . . . . 8 2.4 Direct Numerical Inversion . . . . . . . . . . . . . . . . . . . . . 9 2.4.1 Weak form of forward problem . . . . . . . . . . . . . . . 9 2.4.2 Linearized inverse problem . . . . . . . . . . . . . . . . . 11 3 A Simple Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4 Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 17 4.2 Complex Variables . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.3 Calderón’s Method . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.4 Relation to the Schrödinger Equation . . . . . . . . . . . . . . . 21 4.5 Beals and Coifman’s Method . . . . . . . . . . . . . . . . . . . . 21 4.6 Nachman’s Solution . . . . . . . . . . . . . . . . . . . . . . . . . 25 ii Table of Contents 4.7 Summary of Nachman’s Method . . . . . . . . . . . . . . . . . . 27 4.8 Comparison to Brown and Uhlmann’s Method . . . . . . . . . . 28 5 Approximate Solution . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.1 A Survey of Numerical Implementations . . . . . . . . . . . . . . 30 5.2 The Dirichlet-Neumann Map . . . . . . . . . . . . . . . . . . . . 32 5.2.1 Experimental determination . . . . . . . . . . . . . . . . 32 5.2.2 Radially symmetric conductivities . . . . . . . . . . . . . 33 5.2.3 Single disc . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.3 The Boundary Integral . . . . . . . . . . . . . . . . . . . . . . . 35 5.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.3.2 Fadeev’s Green’s function and its calculation . . . . . . . 35 5.3.3 Single layer operator Sk . . . . . . . . . . . . . . . . . . . 38 5.3.4 The low frequency approximation . . . . . . . . . . . . . 39 5.3.5 Radial symmetry . . . . . . . . . . . . . . . . . . . . . . 40 5.4 The Scattering Transform . . . . . . . . . . . . . . . . . . . . . . 40 5.4.1 Low Frequency Approximation . . . . . . . . . . . . . . . 41 5.5 The ∂-bar Equation . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.5.1 First-order interpretation of µ(x, 0) . . . . . . . . . . . . 42 5.5.2 Solution in terms of convolution . . . . . . . . . . . . . . 43 5.5.3 Centre of radially symmetric problems . . . . . . . . . . 46 5.6 Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.1 Summary of the Approximation . . . . . . . . . . . . . . . . . . 53 6.2 Opportunities for Further Research . . . . . . . . . . . . . . . . 54 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 II Appendices 59 A Schrödinger equation and Laplace’s equation . . . . . . . . . . 60 B Potential for One and Two Discs . . . . . . . . . . . . . . . . . . 61 C Single disc scattering transform, low-frequency approximation 70 iii List of Tables 3.1 Image contrast vs. conductivity . . . . . . . . . . . . . . . . . . . 14 4.1 Comparison of Nachman’s and Brown & Uhlmann’s solution . . 29 5.1 Numerical implementations of Nachman’s or Brown & Uhlmann’s solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 iv List of Figures 2.1 Back-projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.1 Image contrast vs. conductivity . . . . . . . . . . . . . . . . . . . 15 3.2 Image signal strength vs. radius . . . . . . . . . . . . . . . . . . . 16 4.1 Product of complex exponentials forming plane-wave. (a) eikx (b) eikx (c) eikx+ikx . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2 Conductivity cross-section . . . . . . . . . . . . . . . . . . . . . . 22 5.1 Propagation of potential from origin to boundary for radially- symmetric piece-wise constant conductivities. . . . . . . . . . . . 33 5.2 Exponential solutions and Green’s functions. 5.2(a) an exponen- tial solution ψ = eikx and 5.2(b) the corresponding µ = 1. 5.2(c) Fadeev’s Green’s function Gk(x) and 5.2(d) the corresponding gk = e−ikxGk(x). . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.3 ∂-bar solution µ(0, x) for a disc object . . . . . . . . . . . . . . . 45 5.4 Conductivity at origin computed from ∂-bar equation. . . . . . . 47 5.5 Image reconstruction for disc of moderate and high conductivity. 49 5.6 Image reconstruction for annulus of moderate and high conduc- tivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 B.1 A harmonic potential is scattered from a metallic disc. The ap- plied, scattered, and true potentials are shown in B.1(a), B.1(b), and B.1(c) respectively. The metallic disc is an equipotential, and the scattered potential decays. . . . . . . . . . . . . . . . . . . . 64 B.2 Reflected images of two discs. . . . . . . . . . . . . . . . . . . . . 66 B.3 A harmonic potential is scattered from a two discs. The applied, scattered, and true potentials are shown in B.3(a), B.3(b), and B.3(c) respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . 69 v Part I Thesis 1 Chapter 1 Introduction Electrical impedance tomography (EIT) is a method of producing cross-sectional images of an object from measurements of electrical impedance at its boundary. Geophysical applications such as electrical well-logging and galvanic resistivity measurements are examples of this method. There are also some applications in the monitoring of industrial mixing chambers, and human breathing. This field has a long history, with the earliest application dating back to the 1930’s when the Schlumberger brothers, Conrad and Marcel, analyzed electrical sounding data. Their innovations would be used to start Schlumberger Well Services, which has since evolved into the oilfield services company Schlumberger Limited whose market capitalization on August 14th, 2008, exceeds $100B. Far less august are the applications to biomedical and industrial imaging. In the 1980’s Barber and Brown applied a variant of the geophysical technique to monitoring the heart and lungs. Despite its favourable economics and penetrating power, widespread appli- cation has been limited by the generally poor quality of resulting images and the computational difficulty of generating them. To obtain the electrical impedance cross-section, it is necessary to solve a non-linear inverse problem. The stan- dard approach to do this is through iteratively solving a linearized version of the problem, however the solution is ill-posed, effectively requiring the inversion of a matrix with near-zero eigenvalues. Regularization, a method to stabilize small eigenvalues, is necessary. In 1980, Calderón wrote a paper outlining a new strategy for solving the two- dimensional inverse problem. Although he only described a solution to the linear problem, his strategy suggested that a new analytic solution would be possible. It took several years before significant progress was made. Beals and Coifman outlined the solution as part of their research in inverse scattering, summarized in their 1989 review paper[6]. Adrian Nachman proved the inverse scattering result in 1996 and related it to boundary data [33]. Since then there have been numerous papers trying to alleviate the continuity constraints in Nachman’s 2 Chapter 1. Introduction proof. See for example Brown and Uhlmann [9] and Astala and Päivärinta in [3]. 1.1 Overview of Nachman’s Construction Numerically implementing Nachman’s construction is challenging. It involves manipulating equations varying on exponentially wide scales. There are five key steps: (a) Measuring boundary voltages and currents. (b) Solving a boundary-integral for potential solutions with a particular asymp- totic form. (c) Computing a scattering transform from the boundary integral solution. It is related to the fourier transform of the object being imaged. (d) Solving the ∂-bar equation, as described in Section 5.5. This step uses the scattering transform to compute potentials within the domain and away from the boundary. (e) Computing the conductivity within the domain. Measurements of boundary voltages and currents are the starting point for this reconstruction. Nachman’s proof assumes the ability to apply currents to and make measurements from the entire boundary of a closed two-dimensional domain. Although approximations may be made without complete boundary data, all existing work assumes that complete boundary measurements are avail- able. The boundary-integral of the second step is ill-posed and challenging to regularize. Existing numerical work skips this step and instead uses the known solution for a homogeneous medium. However, it will be shown that an improved approximation to this integral can significantly improve image quality. The computation of the scattering transform traditionally involves Taylor series expansion. Instead, it will be shown that the inverse Fourier transform of the scattering transform, i.e., a linearized image, may be used as a basis for image reconstruction. The ∂-bar equation of the fourth step has been well studied in the literature, and a solution to the equation in terms of Fourier transforms is presented. In 3 Chapter 1. Introduction the radially symmetric test disc-shaped and annular regions examined, however, the ∂-bar equation does little to improve image quality. The final step in Nachman’s reconstruction is to compute the conductivity based on solutions to the ∂-bar equation. A new scaling is introduced here that improves upon Nachman’s technique in the context of this approximation (Nachman’s construction is exact without the approximation presented here, but is difficult to implement numerically). 1.2 Organization of Work A survey of of various inversion strategies is presented in Chapter 2. These techniques are not directly related to Nachman’s constructive uniqueness proof, but rather provide breadth from which to judge the results. This is followed by a simple example in Chapter 3 of the imaging of a single disc-shaped object. This example illustrates the ill-posedness of the problem and suggests a numerical scaling for image contrast that is used in the EIT approximation. Chapter 4 is a review of Nachman’s construction. It does not prove the uniqueness theorem, but rather walks through the mechanical steps, including surveys of relevant work by Calderón, Beals, and Coifman. This is essential to understanding Chapter 5, the new approximation. Finally, concluding remarks and opportunities for further research are given in Chapter 6. 4 Chapter 2 Survey of Electrical Impedance Tomography 2.1 Apparent Resistivity The history of electrical impedance tomography EIT is split into two disciplines of application, geophysical and biomedical, as well as two domains of research, applied and theoretical. The technique first developed in the 1930’s for use in geophysics as a means of surveying the material properties of the earth at great depths. Initially, the primary application was studying vertical stratification. R.E. Langer [28] showed how conductivity below the surface could in principle be reconstructed if it were analytic in depth by measuring the potential near a single electrode injecting current at the surface. The potential beneath the plane of the ground would satisfy, (2.1) σ(z) [ ∂2u ∂r2 + 1 r ∂u ∂r + ∂2u ∂z2 ] + dσ(z) dz · ∂u ∂z = 0. Langer used separation of variables to decompose u(r, z) = R(r)Z(z) into R′′ + r−1R′ + λ2R = 0(2.2) {σZ ′}′ − λ2σZ = 0.(2.3) The equation for R is a Bessel equation with solution R(r) = J0(λr), and that for Z is a Sturm-Liouville equation with solutions of exponential form, denoted Z1(x,λ). By requiring ∂u/∂z = 0 at the surface except within a radius a where current I is injected, Langer reasoned that (2.4) u(r, z) = −I 2piaσ(0) ∫ ∞ 0 Z1(z,λ) Z ′1(0,λ) J0(λr) sinλa dλ. 5 Chapter 2. Survey of Electrical Impedance Tomography By inverting the Bessel equation and using a series solution for the Sturm- Liouville equation, he proceeded to solve the inverse problem for analytic con- ductivities. In applied geophysics, however, Langer’s analytical approach was not used because it did not address the instability in the inverse problem, nor handle the discontinuous conductivities typically sought in applications. Instead, the con- cept of apparent resistivity was introduced to bypass the inversion. Stefanesco and Schlumberger [39] expressed the potential around a single current source above a horizontally stratified medium as (2.5) V (r) = I 2pi ρ1 { 1 r + ∫ ∞ 0 K(λ)J0(λr) dλ } , where K(λ) is a reflection kernel depending on the vertical stratification. This equation is much like Langer’s Equation 2.4. The apparent resistivity, ρa at depth h is calculated based on the observed electric field at a distance h from a current source, (2.6) ρa = −(2pi/I)r2∂V/∂r, i.e., the apparent resistivity is the homogeneous resistivity that would generate the observed electric field corresponding to the injected current. One- and two-dimensional plots of apparent resistivity were collected and special cases of horizontal stratification were analyzed by curve-fitting against analytical results [24], [25], [26], [12]. Early results for 1-dimensional stratifi- cation often relied on current injection at a single point, and measurement of the surface potential at various radii as with Langer’s approach. However, there were several experimental difficulties with this approach, and by the 1950’s, many new electrode arrays were developed. For a recent comparison of the ma- jor types see [15]. In the 1980’s Barber and Brown applied the geophysical techniques to med- ical imaging and re-stimulated algorithmic research with a simple 2D back- projection algorithm that generated reasonable images quickly [5], [8]. Barber and Brown considered a circular domain with current source and sink at op- posite ends. For this configuration, equipotentials are circular, even if a strip between two equipotentials has differing conductivity, Figure 2.1. By measur- ing differential voltages along the circumference for a given current pattern, a conductivity value may be assigned to each strip bounded by equipotentials. 6 Chapter 2. Survey of Electrical Impedance Tomography Superimposing the conductivities calculated for every combination of opposite source and sink leads to a blurry image of the interior. This technique is similar to X-ray tomography where densities are projected along the straight photon transmission path. For this reason, Barber named the imaging technique electri- cal impedance tomography. Both the apparent resistivity methods of geophysics !0!1 V Figure 2.1: Back-projection and Barber and Brown’s back-projection algorithm are quick to calculate and stable, however they do not attempt to directly invert the conductivity problem to generate an improved image. By the 1990’s, it became computationally feasible to apply numerical tech- niques such as Tikhonov regularization [40] to compute true resistivity from measurements in two and three dimensions [29]. 2.2 Variational Methods Berryman, Kohn, and McKenney [7] and [27] developed a variational technique that treated ∇ · J = 0 and E = −∇Φ as constraints and minimized the error in J = σE. They did this with a novel functional, I(σ,Φ, J) = ∫ Ω ∣∣∣σ1/2∇Φ+ σ1/2J∣∣∣2 dx(2.7) = 1 2 ∫ Ω σ|∇Φ|2 dx+ 1 2 ∫ Ω σ−1|J |2 dx+ ∫ Ω J ·∇Φ dx.(2.8) 7 Chapter 2. Survey of Electrical Impedance Tomography The first term in the expansion of Equation 2.8 is Dirichlet’s variational principle for ∇ ·(σ∇Φ) = 0. Since this quantity is minimized by reducing σ(x), it provides an upper limit to σ. The second term is Thompson’s dual variational principle and provides a lower limit to σ. The third is calculated from the boundary data and follows from ∇ · J = 0. Although novel, the algorithm’s numerical implementation was simplistic and produced images with significant artifacts. Another variational algorithm for the linear case was developed by Dobson and Santosa in [17]. They sought a conductivity profile with minimum total variation, (2.9) min I(δσ) = ∫ Ω |∇δσ|. This however, excessively penalizes variation near the boundary where resolution is highest. 2.3 Resolution and Optimal Current Patterns Resarch by Cheney and Isaacson [13] defined a notion of distinguishability in EIT. This is important because it allows for solutions to penalize large variation differently depending on position. Resolution scales inversely with depth, so much higher variation can be tolerated near the surface than in the interior. Cheney and Isaacson’s distinguishability criterion is (2.10) δ1/2(σ,σ0, j) = ||v − v0||1/2 ||j||−1/2 , where ||f ||p is the Sobolev norm of order p, (2.11) ||f(x)||p = ||(1 + |k|p)F [f(x)] ||2, and F is the Fourier transform. It can be used to optimize current patterns for simple shapes to maximize SNR. Finding optimal current patterns reduces to calculating the eigenvectors for a system of equations. For the case of a centered cylinder in a circular domain, the optimal patterns are sinusoidal, and larger electrodes that more closely approximate the eigenfunctions lead to improved imaging performance. In practice, it is quite difficult to approach this ideal, and generally impossible for geophysical measurements. 8 Chapter 2. Survey of Electrical Impedance Tomography Ping Hua et al. used an iterative approach to optimizing current patterns in [18], a FEM based iterative Newton-Raphson algorithm, and then extended it for parallel computation [41]. Paulson et. al. [34] added a Robin condition to better model contact impedance. Denoting the boundary of the kth electrode by Γk, the boundary conditions become∫ Γk σ ∂φ ∂n ds = Ik,(2.12) φ+ zσ ∂φ ∂n = Vk,(2.13) where φ is the potential within the domain, σ is the conductivity at the sur- face, Ik is the current through the kth electrode, Vk is the potential at the kth electrode, and z is the contact impedance. The contact impedance is not only a function of electrode geometry, but also a frequency-dependent function of the electrochemistry at the interface. 2.4 Direct Numerical Inversion N. Polydorides andW. Lionheart developed a Matlab toolkit: Electrical Impedance and Diffuse Optical Reconstruction Software (EIDORS) as a testbed for nu- merical algorithms [35]. The software implements a conventional finite element forward solver and some novel regularization techniques for inversion. The weak formulation and Tikhonov regularization are described here. 2.4.1 Weak form of forward problem The Laplace equation with Robin boundary conditions is ∇ · σ∇u = 0(2.14) au+ bσ ∂u ∂ν = g(2.15) Away from the electrodes, Neumann boundary conditions (a=g=0, b=1) ensure that no current crosses the boundary. The electrodes are modeled by a constant 9 Chapter 2. Survey of Electrical Impedance Tomography impedance model that is accurate to within the imaging capabilities of EIT [19]. (2.16) Vl = u+ zlσ ∂u ∂n Il = ∫ Γl ∂u ∂n σ ds 1 < l < L The weak form is obtained by multiplying Equation 2.14 by a test function v and integrating by parts. (2.17) ∫ Γ vσ ∂u ∂n ds− ∫ Ω ∇v · σ dx = 0. (2.18) ∫ Γ v g − au b ds− ∫ Ω ∇v · σ∇u dx = 0. (2.19) L∑ l=1 ∫ Γl v Vl − u zl ds− ∫ Ω ∇v · σ∇u dx = 0. This is a bilinear equation of the form a(u, v) = l(v) where (2.20) a(u, v) = ∫ Ω ∇v · σ∇u dx+ L∑ l=1 ∫ Γl vu zl ds l(v) = L∑ l=1 ∫ Γl v Vl zl ds. If currents rather than voltages are applied to the boundary, Il is related to Vl by (2.21) Il = ∫ Γl Vl − u zl . For there to be a solution, a necessary condition is that (2.22) ∑ l Il = ∫ Γ Vl − u zl = 0, which can be seen by setting v = 1 in Equation 2.19. Physically this means that 10 Chapter 2. Survey of Electrical Impedance Tomography there can be no net charge accumulation within the domain in this steady-state or time-harmonic system. Since a(u, v+ c)− l(v+ c) = a(u, v)− l(v), the system of equations 2.20 and 2.21 have a unique solution to within a constant voltage. 2.4.2 Linearized inverse problem The weak formulation defines a linear relationship between the electrode cur- rents, Il, and the electrode potentials, Vl, parametrized non-linearly by the con- ductivity σ. Consider a set of measured current patterns and voltages, F (σ) : I(k)l → Vl(k), where l is the electrode number and k is the current pattern. In the EIDORS package, adjacent electrodes are driven in sequence around a circular or cylindrical domain. The goal of the inverse problem is to choose σ to minimize ||F (σ) − V ||22. By linearizing about a specific conductivity distribution, σ0, F (σ) = F (σ0) + F ′(σ0)h. The goal becomes to choose h to minimize ||F (σ0) + F ′(σ0)δσ − V ||22. Since the inverse problem is ill-posed, there are infinitely many solutions to this problem and it is necessary to penalize solutions that are unlikely. With Tikhonov regularization, the problem is to minimize (2.23) ||F (σ0) + F ′(σ0)δσ − V ||22 + λ||Hδσ||2, where H is an operator that quantifies the badness of a particular conductivity and λ is a regularization parameter that penalizes the badness. For example, in this thesis (which does not use the EIDORS toolkit), Hδσ ≈ ∆(δσ), so conductivity distributions that are rapidly varying are penalized. Equation 2.23 is discretized in terms of finite elments with piecewise constant conductivity, and the minimum is achieved for h solving (2.24) F ′†F ′δσ + λH†Hδσ = F ′†(V − Fσ0). This method for improving an estimate for conductivity may be iterated, and λ may vary with each iteration. Normally, however, the quality of measured data does not justify more than one iteration, and many algorithms including the EIDORS test cases and Noser [14] stop after the first iteration from a uniform σ0. Most modern inversion algorithms solve the weak forward problem with finite elements or finite differences, and use it to compute the Jacobian F ′ for inversion. 11 Chapter 2. Survey of Electrical Impedance Tomography Tikhonov and other more sophisticated strategies are used for regularization. Nachman’s construction is quite different from this direct inversion method. However, the underlying physical problem is the same, so similar strategies are used to handle the ill-posedness. But first, we examine a simple example. 12 Chapter 3 A Simple Example A single disc object is a useful test case because it easy to solve the forward and inverse problem analytically, and because it shares many features with more sophisticated algorithms for more general conductivities. Consider Laplace’s equation in two dimensions, ∇ · σ∇u = 0, with radially symmetric conductivity. Furthermore, suppose the domain for the problem is the full complex plane, and that the conductivity is σ1 inside a radius of R1 and σ2 outside R1. In each region, the functions { 1, zn, z−n, zn, z−n } satisfy Laplace’s equation and form a basis for solutions. At interface R1, Laplace’s equation reduces to a boundary condition, u|r=R−1 = u|r=R+1(3.1) σ1ur|r=R−1 = σ2ur|r=R+1(3.2) If the interior solution is u1 = 1 · zn the exterior solution is of the form u = azn + bz−n = arneinθ + br−neinθ where, 1 ·Rn1 = a(n)Rn1 + b(n)R−n1 σ1nR n−1 1 = a (n)σ2nR n−1 1 − b(n)σ2nR−n1 By solving this system, the solution in the exterior is a(n) = (1 + σ1σ2 2 ) b(n) = (1− σ1σ2 2 ) R2n1 u = (1 + σ1σ2 2 ) zn + (1− σ1σ2 2 ) R2n1 z −n Consider now the inverse problem: finding σ1 and R1 given measurements of a(1), b(1), a(2), b(2), and a knowledge of σ2. By taking a ratio of the “scattered” 13 Chapter 3. A Simple Example and “incident” components, c(j) = b (j) a(j) we obtain c(1) = ( 1− σ1σ2 1 + σ1σ2 ) R21(3.3) c(2) = ( 1− σ1σ2 1 + σ1σ2 ) R41(3.4) from which R1 = √ c(2) c(1) and the contrast f ≡ ( 1−σ1σ2 1+ σ1 σ2 ) = (c (1))2 c(2) . Note that the a and b measurements represent potentials of the form arneinθ and br−neinθ and are normalized as though the measurements are made at r = 1, i.e. aeinθ, beinθ. So, if one is measuring an object with R1 ' 1, then one must resolve very small scattered potentials b(n) = ( 1−σ1σ2 2 ) R2n1 . These terms decay geometrically with the size of the object. Furthermore note that the contrast term f saturates. If the object matches the background, σ2 = σ1 then f = 0; if the object is a perfect conductor, σ1 = ∞ then f = −1; if the object is a perfect insulator, σ1 = 0 then f = +1. This is illustrated in Table 3.1 and Figure 3.1. Table 3.1: Image contrast vs. conductivity Conductivity Ratio Contrast 0 1 1/7 3/4 1/3 1/2 3/5 1/4 1 0 5/3 −1/4 3 −1/2 7 −3/4 ∞ −1 In a physical measurement, one would not directly measure a and b, nor would one consider the complex plane as a domain. A (slightly) more realistic problem would be to consider a unit disc as the domain (an idealized body), to apply a current pattern j = j(n)einθ, and then measure a voltage pattern 14 Chapter 3. A Simple Example 10 !2 10 !1 10 0 10 1 10 2 !1 !0.5 0 0.5 1 Image Contrast vs. Conductivity Ratio ! 1 / ! 2 (! 1 ! ! 2 ) / (! 1 + ! 2 ) Figure 3.1: Image contrast vs. conductivity v(n)einθ. In this case, v(n) = a(n) + b(n) j(n) = σ2|n| ( a(n) − b(n) ) . Hence, a(n) = 1 2 [ v(n) + j(n) σ2|n| ] b(n) = 1 2 [ v(n) − j (n) σ2|n| ] If σ1 were equal to σ2 then one would expect v = j (n) σ2|n| , a = j(n) σ2|n| , and b = 0. The inverse problem is solved by examining the difference between the measured potentials and the ideal constant conductivity potentials, c(n) = v − j(n)σ2|n| v + j (n) σ2|n| , so accurate knowledge of σ2 is critical. In a more realistic scenario many more measurements would be made including higher frequency terms, e.g. with closely 15 Chapter 3. A Simple Example 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Signal Strength vs. Object Radius R 1 R 12 , R 14 Figure 3.2: Image signal strength vs. radius spaced electrodes, to resolve σ2 and some assumptions would be made about the regularity of the conductivity near the surface. To see that determining σ and R is an ill-posed problem, note that their calculation requires a precise determinatino of a and b and that the b terms vary as R2n. For objects far from the boundary, these values are very small, as shown in Figure 3.2. 16 Chapter 4 Analytical Solution Moving from the simple test case to the general problem, complexity increases significantly. In pure mathematics, the 2-dimensional EIT problem is considered solved with a constructive proof, however there has not been a single complete numerical implementation of this construction. Nachman’s paper on this topic [33] is quite technical, so a a re-derivation is given here. Note that this chapter does not mathematically prove Nachman’s result, but rather motivates the im- portant steps in a manner useful for understanding the approximations in the rest of the thesis. For the remainder of the thesis, the EIT algorithm is specialized to the 2-dimensional problem with circular domain and unit conductivity near the boundary. For a discussion of why this is reasonable, and the relationship be- tween this special case and more general problems, see [33], [38], and [20]. 4.1 Mathematical Formulation Structures with different conductivities modify the surrounding potential in a manner that satisfies Laplace’s equation for inhomogeneous conductivity. (4.1) ∇ · σ∇u = 0, in Ω Impedance measurements are made on the boundary, often by forcing a current through the medium. The electrical impedance tomography algorithm assumes that the impedance at the surface is known for all possible meausrement con- ditions. One way of representing this is as follows: for any given surface voltage f , the current g flowing through the surface is known. This is equivalent to knowing the Dirichelet-Neumann map Λσ, where( σ ∂u ∂n ) = Λσ(f), u|∂Ω = f. (4.2) 17 Chapter 4. Analytical Solution In the study of partial differential equations, one traditionally attempts to solve the forward problem for Equation 4.1. Given a conductivity and a set of boundary conditions, one finds the potential throughout the domain (4.3) σ, f → u in Ω. This is well posed, and there are many efficient algorithms for its solution. The problem of electrical impedance tomography is the inverse, given a set of boundary measurements, to find the conductivity of the interior of the domain, (4.4) Λσ → σ. This problem is non-linear and ill-posed. 4.2 Complex Variables As a notational convenience, two operators for complex variables are used, ∂ and ∂. ∂ = ∂ ∂x = 1 2 ( ∂ ∂x1 − i ∂ ∂x2 ) ∂ = ∂ ∂x = 1 2 ( ∂ ∂x1 + i ∂ ∂x2 ) The reason for these definitions is that for holomorphic functions such as xn where x is complex, (4.5) ∂xn = nxn−1, n ≥ 1 and for antiholomorphic functions (4.6) ∂xn = 0, n ≥ 0. As a result, ∂ commutes with antiholomorphic functions, e.g., (4.7) ∂xnxn = xn∂xn = nxnxn−1 The following elementary results are written for reference. Note that the 18 Chapter 4. Analytical Solution same expressions are valid if their complex conjugate is taken. The Laplacian in terms of these complex derivatives is (4.8) ∆ = 4∂∂ = 4∂∂. The following formula contains the Green’s function for the Laplacian, (4.9) ∆ [ 1 2pi log |x| ] = 4∂∂ [ 1 2pi log(x) ] = δ(x). and the next the Green’s function for the ∂ operator, (4.10) ∂ ( 1 pix ) = δ(x). Additionally, (4.11) ∂ log |x| = 1 2x , where the factor of 1/2 appears because log |x| contains a holomorphic and anti-holomorphic component. 4.3 Calderón’s Method Caldéron showed how to solve for conductivities that are nearly constant [11]. The divergence theorem relates integrals of the spatially varying conductivity to voltage and current measurements on the boundary. 0 = ∫ Ω v∇ · σ∇u dx = ∫ ∂Ω vσ ∂u ∂n ds− ∫ Ω ∇v · σ∇u dx,(4.12) ∫ Ω ∇v · σ∇u dx = ∫ ∂Ω vΛσ(f) ds.(4.13) Consider a conductivity σ = 1 + δσ where the perturbation is strictly within the domain Ω. Let u = eikx + δu and v = eikx, with δu satisfying ∇ · σ∇(eikx + δu) = 0, δu = 0 on the boundary. This choice of u and v makes their product approximately a plane wave, eikx+ikx = ei(2k1,−2k2)·(x1,x2), as 19 Chapter 4. Analytical Solution (a) (b) (c) Figure 4.1: Product of complex exponentials forming plane-wave. (a) eikx (b) eikx (c) eikx+ikx shown in Figure 4.1. Calderon’s Equation 4.13 becomes∫ ∂Ω eikxΛσ(eikx) ds = ∫ Ω δσ∇eikx ·∇ (eikx + δu) dx(4.14) ≈ 2|k|2 ∫ Ω δσeikx+ikx dx(4.15) = 2|k|2 ∫ Ω (δσ)e−i(−2k1,2k2)·x dx(4.16) Setting ξ = −2k yields (4.17) ∫ ∂Ω ei ξx 2 Λσ(ei ξx 2 ) ds = 1 2 |ξ|2 ∫ Ω (δσ)e−iξ·x dx. This shows that the Fourier transform of δσ may be approximated by mea- surements on the boundary. Since a k = 0 measurement (constant voltage) does not result in current, it yields no information. To properly resolve conductivity in a domain, it is necessary to know the conductivity at the boundary. Alternate formulation in terms of moments Allers and Santosa implemented Calderon’s algorithm by reducing the inverse problem to a moment problem and expanding the conductivity in terms of 20 Chapter 4. Analytical Solution Zernike polynomials [2]. These polynomials are related to the Jacobi polynomials by (4.18) Qm,k(r) = P (0,m) k (2m− r). Instead of choosing potentials of the form eikx and eikx, the authors choose rkeikθ/k and rle−ilθ/l. This corresponds to measuring a Fourier series of the boundary potential perturbation δfk(θ) due to a current pattern eikθ. The re- sulting moment problem is (4.19) ∫ 1 0 ∫ 2pi 0 δσ(r, θ)rk+l−1ei(k−l)θ dθ dr = −1 2 ∫ 2pi 0 δfk(θ)e−ilθ dθ. The Zernike polynomials are chosen as a basis because they are orthogonal with weight function rm. The zeros of these polynomials accumulate near r = 1, allowing for higher resolution near the boundary. 4.4 Relation to the Schrödinger Equation Nachman proved the uniqueness of the 2-dimensional inverse problem by exam- ining the properties of an equivalent Schrödinger equation. If the conductivity is twice differentiable, then the Laplace equation, (4.20) ∇ · (σ∇v) = 0, may be transformed into an equivalent Schrödinger equation, (4.21) (−∆+ q)u = 0, where u = σ1/2v and q = ∆σ1/2/σ1/2. The details of this calculation are in Appendix A. Figure 4.2 sketches a cross- section of a discontinuous conductivity, as well as ∆σ1/2 which appears in the product qu. 4.5 Beals and Coifman’s Method Beals and Coifman outlined a strategy for solving the inverse problem for the Schrödinger equation [6]. Although they did not prove uniqueness for the 2- 21 Chapter 4. Analytical Solution !! 1 f(x)=!! f´(x) f˝(x) Figure 4.2: Conductivity cross-section dimensional case, their general ∂-bar method was used by Nachman in his proof. Beals and Coifman, like Calderón, looked for exponentially growing solu- tions. (4.22) u = µ(x, k)eikx ≡ ψk(x), x, k ∈ C, where µ(x, k) ∼ 1 as x→∞. By substituting 4.22 in 4.21, Schrödinger’s equation for µ becomes ∆u− qu = 4∂∂µeikx − qµeikx = 0[−4∂∂ − 4ik∂]µ(x, k) = −qµ(x, k) This equation may be written (4.23) Lkµ(x, k) = −qµ(x, k), where Lk = (−∆− 4ik∂). The Green’s function for the last equation is the solution to (4.24) Lkgk(x) = δ(x), 22 Chapter 4. Analytical Solution and Equation 4.23 may be rewritten (4.25) µ(·, k) = 1− gk ∗ qµ. The free-space Green’s function gk for the linear operator of Equation 4.24 may be calculated by using Fourier transforms. By defining the 2-dimensional transform pair f̂ = ∫ R2 e−ix·ξf(x) dx f = 1 4pi2 ∫ R2 eix·ξ f̂(ξ) dξ, applying it to Equation 4.24,∫ R2 e−ix·ξ [−4∂∂ − 4ik∂] gk(x) dx = 1 = [ |ξ|2|+ 4ik−iξ1 + ξ2 2 ] ∫ R2 e−ix·ξgk(x) dx = 1, and inverting the Fourier transform, one arrives at the expression gk = 1 4pi2 ∫ R2 eix·ξ 1 |ξ|2 + 2k(ξ1 + iξ2) dξ = 1 4pi2 ∫ R2 eix·ξ 1 ξ(ξ + 2k) dξ. Beals and Coifman looked at the properties of the ∂/∂k derivative of Equa- tion 4.25. By rewriting the equation as (4.26) µ(·, k) = [I + gk ∗ (q•)]−1 1, and noting that q is only a function of x, not ξ, one finds from Equation 4.25 that ∂kµ = −gk ∗ (q∂kµ)− [ ∂k(gk∗) ] (qµ),(4.27) ∂kµ = [I + gk ∗ (q·)]−1 [ ∂k(gk∗) ] qµ.(4.28) After applying Equation 4.10, the term associated with ∂gk∗ simplifies consid- 23 Chapter 4. Analytical Solution erably, ∂ ∂k (gk ∗ f) = ∂ ∂k [ 1 4pi2 ∫ R2 eix·ξ f̂(ξ) ξ(ξ + 2k) dξ ] = 1 4pi2 ∫ R2 eix·ξ f̂(ξ) ξ ∂ ∂k 1 (ξ + 2k) dξ = 1 4pi2 ∫ R2 eix·ξ f̂(ξ) ξ 2piδ(2k + ξ) dξ = ei(−2k1x1+2k2x2) 4pik f̂(−2k1, 2k2), As a notational convenience, the exponential term is rewritten as ek(x) = e2k1x1−2k2x2 = ei(kx+kx), and the generalized Fourier transform of qµ(x, k) is defined to be the scattering transform, t(k) = [qµ(·, k)] (̂−2k1, 2k2) = ∫ Ω eikxqµ(x, k)eikxd2x. With these definitions, (4.29) ∂k (gk ∗ (qµ)) = e−k(x)t(k)4pik . Another useful observation is that Lke−kµ = (−∆− 4ik∂k) e−k(x)µ = (−4∂∂ − 4ik∂k) e−ikx−ikxµ = e−ikx−ikx (−∆+4 ik̄∂k)µ = e−kLkµ. This property extends to gk∗ as well, since gk∗ is the inverse of Lk, e−k(µ− 1) = gk ∗ Lk(e−kµ) = gk ∗ e−kLkµ, = e−kgk ∗ Lkµ, 24 Chapter 4. Analytical Solution which implies that gk ∗ e−kµ = e−kgk ∗ µ. Applying this result to Equation 4.26 yields (4.30) [I + gk ∗ (q•)]−1 e−k = e−k [ I + gk ∗ (q•) ]−1 1 = e−kµ(x, k). Combining 4.28, 4.29, and 4.30 gives the result of Beals and Coifman’s ap- proach to this problem, (4.31) ∂kµ(x, k) = e−k(x)t(k) 4pik̄ µ(x, k). This so called ∂-bar equation may be recast as an integral-equation by noting that the inverse of ∂k is 1/pik. This gives an analogue in k to the equation for µ in x, µ(x, k) = [I + gk ∗ (q•)]−1 1 µ(x, k) = [ I − 1 pik ∗ e−k(x)t(k) 4pik̄ •̄ ]−1 1. 4.6 Nachman’s Solution Beals and Coifman only suggested a strategy for solving the Schrödinger equa- tion ∆u−qu = 0. Nachman developed a constructive proof requiring q ∈ Lp, 1 < p < 2, µ(·, k)− 1 ∈ Lp̃ ∩L∞, p̃ = 1/p− 1/2. This essentially requires two deriva- tives in conductivity and that µ be bounded and decay to 1 as x → ∞. Nach- man applied the general method to the EIT problem and proved that boundary measurements uniquely determine conductivity. In terms of the reconstruction algorithm, Nachman made two important observations. (a) The scattering transform, t(k) is uniquely determined from boundary mea- surements. (b) The conductivity may be recovered from µ(x, 0). Let u be a solution to the ∆u − qu = 0 on Ω, w be a harmonic function, and u0 be a harmonic function that is equal to u on ∂Ω. Then applying the 25 Chapter 4. Analytical Solution divergence theorem, ∫ Ω wqu = ∫ Ω w∆u = ∫ ∂Ω w ∂u ∂n − ∫ Ω ∇w ·∇u = ∫ ∂Ω w ∂u ∂n − ∫ ∂Ω u ∂w ∂n + 0 = ∫ ∂Ω w ∂u ∂n − ∫ ∂Ω w ∂u0 ∂n = ∫ ∂Ω w (Λq − Λ0)u, Applying this to the definition of t(k) yields, t(k) = ∫ Ω eikxqµ(x, k)eikx dx = ∫ ∂Ω eikx (Λq − Λ0)µ(x, k)eikx dσ. (4.32) A similar analysis may be used to find µ along the boundary in terms of a single-layer potential. The function Gk(x) = eikxgk(x) solves −∆u = δ(x) and is harmonic everywhere except the origin. µ(x, k)eikx = eikx − ∫ Ω Gk(x− y)q(y)µ(y, k)eiky dy = eikx − ∫ ∂Ω Gk(x− y) (Λq − Λ0)µ(y, k)eiky dσ(y). Defining the operator Sk by (4.33) (Skφ) (x) = ∫ ∂Ω Gk(x− y)φ(y) dσ(y), the equation for µ along the boundary becomes (4.34) µ(x, k)|∂Ω = 1− e−ikxSk (Λq − Λ0) eikxµ(x, k). The scattering transform t(k) may be recovered from boundary measurements by equations 4.32 and 4.34. 26 Chapter 4. Analytical Solution Nachman’s second observation is that µ(x, 0) = σ1/2. This is easily verified by substituting σ1/2 in Schrödinger’s equation. (4.35) (−∆+ q)µ(x, 0) · 1 = (−∆+ ∆σ 1/2 σ1/2 )µ(x, 0) = 0. Nachman’s method has been criticized because of its reliance on twice- differentiable conductivities, however this is a limitation of the theoretical meth- ods involved in the proof. If one considers a piece-wise continuous conductivity, then q = ∆σ1/2/σ1/2 will be undefined, however the product qu = ∆σ1/2v does have a valid interpretation. Furthermore, one may conjecture from the linear problem that since t(k) is related to the Fourier transform of q, t(k)/|k|2 is related to the Fourier transform of σ and decays as k → ∞, for conductivities with jump discontinuities. 4.7 Summary of Nachman’s Method The full conductivity reconstruction procedure is as follows, (a) Measure the Dirichlet-Neumann map Λq (4.36) Λq : u|∂Ω → ∂u ∂n ∣∣∣∣ ∂Ω . (b) Solve for µ(x, k) on ∂Ω (4.37) µ(x, k)|∂Ω = 1− e−ikxSk (Λq − Λ0) eikxµ(x, k). (c) Calculate the scattering transform t(k) (4.38) t(k) = ∫ ∂Ω eikx (Λq − Λ0)µ(x, k)eikx dσ. (d) Solve for µ(x, k) within domain Ω (4.39) µ(x, k) = 1 + 1 pik ∗ e−k(x)t(k)µ(x, k) 4pik . (e) Calculate the conductivity within the domain Ω (4.40) σ(x) = µ(x, 0)2. 27 Chapter 4. Analytical Solution 4.8 Comparison to Brown and Uhlmann’s Method After Nachman solved the inverse conductivity problem in terms of an equivalent second-order differential equation [33], Brown and Uhlmann solved the problem in terms of coupled first-order differential equations [10]. This proof required less regularity in the conductivity. A brief comparison of the two formulations is given in Table 4.1. The reconstruction formulas for Brown and Uhlmann’s solution are from Knudsen [20]. 28 Chapter 4. Analytical Solution T ab le 4. 1: C om pa ri so n of N ac hm an ’s an d B ro w n & U hl m an n’ s so lu ti on E qu at io n N ∆ v − qv = 0 B U ( ∂ 0 0 ∂ )( v w ) − ( 0 q q 0)( v w ) V ar ia bl es N v = σ 1 / 2 u q = ∆ σ 1 / 2 σ 1 / 2 B U ( v w) = σ 1 / 2 ( ∂u ∂u) q = −σ 1 / 2 ∂ σ 1 / 2 A sy m pt ot ic fo rm N ψ (x ,k ) = µ (x ,k )e ik x µ → 1 as |x| → ∞ B U ψ (x ,k ) = m (x ,k )( ei k x 0 0 e− ik x ) m → I, as |x| → ∞ G re en ’s fu nc ti on N G k (x ) = − 1 2 pi .[ E i( ik x )] B U G k (x ) = 1 pi e − i k x x B ou nd ar y in te gr al N µ (x ,k )| ∂ Ω = 1 − e− ik x S k (Λ q − Λ 0 )e ik x µ (x ,k ) B U 1 + iS k 0 0 1 − iS k i( H γ − 1) ν −i (H γ + 1) ν ( ψ 1 1 (x ,k ) ψ 2 1 (x ,k )) = 2ei k x 0 0 Sc at te ri ng tr an sf or m N t( k ) = ∫ ∂Ωe ik x (Λ q − Λ 0 )µ (x ,k )e ik x dσ B U t( k ) = i 2 pi ∫ ∂Ω( 0 e− ix k ν (x )Ψ 1 2 (x ,k ) −e ix k ν (x )Ψ 2 1 (x ,k ) 0 ) dσ ∂ -b ar eq ua ti on N ∂ k µ (x ,k ) = e − k (x )t (k ) 4 pi k̄ µ (x ,k ) B U ∂ k m (x ,k ) = m (x ,k )( e k (x ) 0 0 e − k (x )) t( k ) C on du ct iv it y N σ 1 / 2 = µ (x ,0 ) B U σ 1 / 2 = .[ m̃ + (x ,0 )] 29 Chapter 5 Approximate Solution 5.1 A Survey of Numerical Implementations There have been several attempts to numerically implement either Nachman’s (N) or Brown and Uhlmann’s (BU) solution to the inverse conductivity problem. Each implementation has had its successes and limitations as outlined in Ta- ble 5.1. All have assumed an idealized circular domain with unit conductivity at the boundary and sinusoidal stimuli. Although this appears to be a limitation, it is shown in [33] that this deliberate simplification does not materially alter the problem.1 It is possible to transform a problem with non-unit conductivity near the boundary to an equivalent one with unit conductivity. Furthermore, Nachman’s construction is defined for more general closed domains. Some of the implementations have been specialized to radially symmetric conductivities, while others are two dimensional. This too is a simplification which does not materially affect the algorithms, but does limit an examination of artifacts in non-radially symmetric problems. In fact, a radial problem only simplifies the calculation of the scattering transform t(k), as the solution to the ∂−bar equation is inherently two-dimensional, even for radial problems. Early implementations focussed on C∞ smooth conductivities for two rea- sons: Nachman’s proof required two derivatives, and early approximations to t(k) involved a truncated Taylor expansion about k = 0 and as a result suffered from a small radius of convergence. Knudsen’s early work, [20], [21], [30], used C1 smooth surfaces because of a restriction in Brown and Uhlmann’s proof. Since then, however, discontinuous surfaces [22] have been attempted. Existing algorithms all suffer from a limitation on the effective range of con- ductivity. If an image has very high contrast, modeling errors will introduce large artifacts. Current approximations allow for a range of approximately 1/2 to 2 times the boundary conductivity. The reason for this limitation is not due 1In the post-processing step of this work, Section 5.6, the boundary conductivity does indeed have an effect. 30 Chapter 5. Approximate Solution to Nachman’s or Brown & Uhlmann’s reconstruction procedure, but rather to a coarse approximation in the boundary integral equation 4.37. All the exist- ing algorithms use µ(x, k)|∂Ω = 1. The reason for this approximation is that there are significant numerical difficulties in the solution of the boundary in- tegral. This work improves upon existing approximations by solving 4.37 with S0 instead of Sk, i.e. by making a low-frequency approximation. As a result, the algorithm tolerates a much wider range of conductivities. For disc shaped objects, conductivities between 0 and ∞ are perfectly resolved (limited only by measurement accuracy and the ill-posedness). Early implementations approximated the scattering transform with a Taylor series, (5.1) texp(k) = ∫ ∂Ω eikx (Λq − Λ0) eikx ds = ∑ tmnk m kn, |k| < kcut For radial problems, this method has the advantage that tmn are directly related to boundary measurements made with sinusoidal stimuli, but suffer from the disadvantage of poor Taylor series convergence. Knudsen [20] improved upon this for the special case of a disc scatterer by expanding in terms of a truncated Bessel series, (5.2) texp(k) = −4piα|k|r ∞∑ m=0 (αr)mJ1(2rm+1|k|), where α = σ−1σ+1 and r is the radius of the disc. Note that this is only for testing the algorithm against a known test case, and is not applied to more general conductivities. The Taylor/Bessel approximation is improved in this thesis. Instead of calculating t(k), an equivalent linearized conductivity c(x) is computed such that the Fourier transform of c(x) evaluated at (−2k1, 2k2) gives t(k). The distribution c(x) is found from the solution to an ill-posed linear problem with Tikhonov regularization. The ∂-bar equation 4.39 has been solved in a variety of ways in the literature, including a full matrix solution, 1- and 2-grid methods, and a simple FFT. This work independently derives an FFT-based solution. Finally, this work differs from existing algorithms in that it adds a post- reconstruction scaling that provides an exact solution to disc-shaped objects. 31 Chapter 5. Approximate Solution Table 5.1: Numerical implementations of Nachman’s or Brown & Uhlmann’s solution Paper Sol. Dim. Smooth. Range Scat. Trunc. ∂-bar Siltanen 2000 [37] N Rad C∞ 1-4 q,Λeikx Tayl Mat Mueller 2000 [31] N 2D C∞ 0.5-1.2 q - Mat Knudsen 2002 [20] BU Rad C1.1 0.6-1.3 Λeikx Tayl FFT Knudsen 2003 [21] BU Rad C1.1 0.6-1.3 Λeikx Bess FFT Knudsen 2004 [23] BU 2D C1.1 0.5-2.5 q Bess 2Gr. Mueller 2004 [30] N 2D Exp. 365-426 Λeikx Bess 2Gr Knudsen 2006 [22] N Rad disc. 1-8 a Λeikx Bess FFT Pereira 2008 N Rad C0 0-∞ Λψ(0) Comp FFT aImage quality is poor with higher conductivities. 5.2 The Dirichlet-Neumann Map 5.2.1 Experimental determination The first step in Nachman’s reconstruction method is to measure the Dirichlet- Neumann map, Equation 4.36 which is reproduced here for convenience, Λq : u|∂Ω → ∂u ∂n ∣∣∣∣ ∂Ω . For the case of σ = 1 near the boundary, the potential in Laplace’s equation, uL, is equal to the potential in Schrödinger’s equation uS = σ1/2uL = uL, and so ∂uSdn = ∂uL dn and Λq = Λσ. The computation for t(k) hinges on the knowledge of (Λq − Λ0). Experimen- tally, this is equivalent to knowledge of the difference between measured currents and that expected for a homogeneous domain given a particular voltage stimu- lus. Boundary voltages and currents are approximated by sinusoids in order to simplify the solution. For this thesis, it is assumed that sinusoidal currents may be applied to the boundary and that the Fourier series of voltages be measured. More general boundary conditions with realistic electrodes are the subject of further research and require an additional approximation, particularly with their effect on the boundary integral equation. 32 Chapter 5. Approximate Solution 5.2.2 Radially symmetric conductivities The Dirichlet-Neumann map for radially symmetric, piecewise constant conduc- tivities may be found analytically, and therefore serves as a useful test-case. In each radial segment, the conductivity may be expanded as u(x) = ∞∑ n=1 unx n + vnxn. To find the map, one searches for solutions that (a) satisfy u(x) = einθ|∂Ω, (b) satisfy Laplace’s equation at each radial discontinuity, (c) do not have a pole at the origin. [ 1 0 ] [ a2 b2 ] [ aN bN ] . . .σ1σ2. . .σN Figure 5.1: Propagation of potential from origin to boundary for radially- symmetric piece-wise constant conductivities. Denote the radii of the segments r0, . . . , ri, . . . rN where r0 = 0 and rN = 1, and denote the conductivity of each segment by σ1, . . . ,σi, . . . ,σN . For the innermost segment, there is no pole. Define the potential to be u(r, θ) = 1 · 33 Chapter 5. Approximate Solution r|n|einθ. In each subsequent segment, the potential is expressible as ui(r, θ) = a(n)i r |n|einθ + b(n)i r −|n|einθ. Figure 5.1 illustrates the labeling of conductivity and potential in each region. The solution to Laplace’s equation may be propagated from the origin by enforcing the boundary conditions at the edges of the segment. Both u and σ ∂u∂r are continuous, (5.3) [ Rni R −n i σi+1Rni −σi+1R−ni ][ ai+1 bi+1 ] = [ Rni R −n i σiRni −σiR−ni ][ ai bi ] . This may be solved for a matrix which propagates the solution from one region to the next, (5.4) [ ai+1 bi+1 ] = ( I + σ2 σ3 − 1 2 [ 1 R−2n −R2n 1 ])[ ai bi ] . The potential and normal derivative at RN = 1 are u|∂Ω = aN + bN(5.5) ∂u ∂n ∣∣∣∣ ∂Ω = |n| (aN − bN ) ,(5.6) so the the Dirichlet-Neuman map for a radial piecewise constant conductivity may be defined with a Fourier basis, (5.7) Λqeinθ = |n| ( aN − bN aN + bN ) , which differs from the unit conductivity map by (5.8) (Λq − Λ0) einθ = − 2|n|bN aN + bN . 5.2.3 Single disc For the special case of a single disc of conductivity σ and radius R embedded in a larger disc of unit conductivity and unit radius, (5.9) (Λq − Λ0) einθ = |n| ( 1− 1− cR 2|n| 1 + cR2|n| ) , where c = 1−σ1+σ . 34 Chapter 5. Approximate Solution 5.3 The Boundary Integral 5.3.1 Overview The next step in the reconstruction procedure is to solve the boundary integral equation 4.37, µ(x, k)|∂Ω = 1− e−ikxSk (Λq − Λ0) eikxµ(x, k). The single-layer operator Sk is a convolution with Fadeev’s Green’s function Gk. Prior to the boundary integral’s solution, it is useful to rewrite the Gk in terms of simpler functions that are easier to compute. 5.3.2 Fadeev’s Green’s function and its calculation Fadeev’s Green’s function is real and harmonic except at the origin, and is defined by (5.10) Gk(x) = eikx (2pi)2 ∫ R2 eix·ξ ξ(ξ + 2k) d2ξ. It may be expressed in terms of a more common transcendental functions by first evaluating its derivatives, and comparing its expansion to that of the complex Exponential Integral function [1]. The derivative with respect to k is ∂Gk ∂k = eikx (2pi)2 ∫ R2 ∂ ∂k ei ξx+ξx 2 2ξ ( ξ 2 + k ) d2ξ =eikx ∫ R2 ei ξx+xξ 2 8piξ δ ( k + ξ 2 ) d2ξ =eikx ∫ R2 ei ξx+xξ 2 2piξ δ ( ξ + 2k ) d2ξ =− e ikx 4pik (5.11) 35 Chapter 5. Approximate Solution Similarly, the derivatives with respect to k,x, and x may be evaluated. ∂Gk ∂k = −e ikx 4pik ∂Gk ∂x = −e −ikx 4pix ∂Gk ∂x = −e ikx 4pix (5.12) By integrating the series expansion of the real and harmonic parts, one gets a series expansion that is convergent everywhere but the origin, Gk = ∂−1k ∂kGk + ∂ −1 k ∂kGk + f(x) =− ∂−1k ∞∑ n=0 (ikx)n 4pin!k − ∂−1 k ∞∑ n=0 (ikx)n 4pin!k + f(x) =− log |x| 2pi − ∞∑ n=1 (ikx)n + (ikx)n 4pinn! + f(x). (5.13) The complex exponential integral function has the following expansion. (5.14) Ei(z) = γ + log z + ∞∑ n=1 zn nn! , and has the asymptotic behaviour Ei(z)e−z < −1/z. Fadeev’s Green’s function decays as x → ∞ and is related to the real part of the complex exponential integral function by Gk(x) = − 12pi. [Ei(ikx)] = − 1 2pi [ γ + log |kx|+ ∞∑ n=1 (ikx)n + (ikx)n 2nn! ] . (5.15) A plot of Gk and the exponentially weighted gk = e−ikxGk(x) is in Figure 5.3.2. This relationship between Fadeev’s Green’s function and Ei appears to be new. 36 Chapter 5. Approximate Solution (a) (b) (c) (d) Figure 5.2: Exponential solutions and Green’s functions. 5.2(a) an exponential solution ψ = eikx and 5.2(b) the corresponding µ = 1. 5.2(c) Fadeev’s Green’s function Gk(x) and 5.2(d) the corresponding gk = e−ikxGk(x). 37 Chapter 5. Approximate Solution 5.3.3 Single layer operator Sk The single layer operator Sk of Equation 4.33 may be evaluated on the perimeter of the unit disc in terms of a Fourier expansion. (Skφ) (x) = ∫ ∂Ω Gk(x− y)φ(y) dσ(y) = ∫ ∂Ω ∞∑ m=−∞ ∞∑ n=0 − 1 2pi ×[ γ + log |k(x− y)|+ ∞∑ n=1 (ik(x− y))n + (ik(x− y))n 2nn! ] φ̂mym 2pi dσ(y), (5.16) where φ̂m = ∫ ∂Ω e −imθφ(θ) dθ. This integral may be evaluated termwise by letting x = einθ and y = einθ ′ and using (5.17) ∫ 2pi 0 log |ik(eiθ − eiθ′)|einθ′ dθ′ = −e inθ 2|n| ∫ 2pi 0 ( ik(eiθ − eiθ′) )n 2nn! e−imθ ′ dθ′ = (−1)m+1(ik)n 2nm! (n−m)!e i(n−m)θ n ≥ m,(5.18) as well as the complex conjugate equation. In principle, the boundary integral equation may be solved in terms of a Fourier expansion for ψk(x) = µ(x, k)eikx on the boundary of the unit disc, (5.19) ψ(n+1)k (x) = e ikx − Sk (Λq − Λ0)ψ(n)k (x), ψ(0)k = eikx, where ψk on the boundary is expanded as (5.20) ψk(x) = n=∞∑ n=−∞ m1=∞∑ m1=0 m2=∞∑ m2=0 ψn,m1,m2e inθ(ik)m1(ik)m2 . Unfortunately, there were stability issues in implementing it numerically. Solving for ψ(x, k) or µ(x, k) on the boundary is similar to solving for the conductivity everywhere. The equation needs regularization, which is most easily provided by a Fourier transformation of µ(x, k) into µ(x, y). In terms of a spatial variable y, 38 Chapter 5. Approximate Solution µ(x, y) is closely related to the gradient of σ(y−x), along with a set of decaying images (see 5.5.2). Unless a constraint on the support of σ(y − x) is provided, the boundary integral may not converge. 5.3.4 The low frequency approximation In comparison to the boundary integral with Fadeev’s kernel above, the log- arithmic term in the Green’s function expansion allows for a straight-forward approximation to ψ(x, k). Define the operator S0 by (5.21) (S0φ) (x) = − 12pi ∫ ∂Ω log |x− y|φ(y) dσ(y) This function may be used in place of Sk to form the approximation. This approximation has several useful features. The operator S0 is indepen- dent of k, and when substituted in the linear boundary integral equation 4.37, the entire equation becomes independent of k except for the inhomogeneous term, (5.22) ψ(x, k)|∂Ω = eikx − S0 (Λq − Λ0)ψ(x, k). The equation may be solved by Taylor-expanding the inhomogeneous component and solving term-wise, (5.23) ψ(n)(x, k)|∂Ω = (ikx) n n! − S0 (Λq − Λ0)ψ(n)(x, k). The solution for ψ0(x, k) is just ∑∞ n=0 ψ (n)(x, k). If there is no object to be imaged, (Λq − Λ0) = 0, and ψ(n) = (ikx) n n! . Even more generally, for any applied potential ψf , one may solve for a func- tion f satisfying (5.24) ψf (x)|∂Ω = f(x)− S0 (Λq − Λ0)ψf (x), i.e. f = [1 + S0 (Λq − Λ0)]ψf . If the function f on the boundary is extended to an entire function in the plane, vb(x), it has an interpretation as an unperturbed background potential. Extending the conductivity distribution to be 1 outside of the domain, the Schrödinger potential ψ(x) is equal to the Laplace potential v(x). It may be thought of as the sum of a background potential and scattered term due to 39 Chapter 5. Approximate Solution objects within the domain, ψ(x) = vb+vs. In the low-frequency approximation, the logarithmic kernel results in vs being harmonic and decaying as x → ∞. Without the low-frequency approximation, vseikx would decay as x → ∞. See Appendix B for a calculation of vb and vs for scattering from one and two discs. 5.3.5 Radial symmetry For radially symmetric problems, the integral equation takes on a particularly simple form as there is no coupling between any of the terms, (5.25) ψ(n)k = (ikx)n n! − ( − 1 2|n| ) (Λ− Λ0)ψ(n)k Denoting (Λq − Λ0)xn = Lnxn, ψ(n)k = (ikx)n n! ( 1− Ln2|n| )(5.26) Note that xn in these equations is equivalent to einθ at r = 1, and does not directly relate to the analytic continuation of the functions to the exterior or interior of the domain. 5.4 The Scattering Transform After obtaining Λ−Λ0 and µ(x, k)|∂Ω one may formally compute the scattering transform t(k) through Equation 4.38. t(k) = ∫ ∂Ω eikx (Λq − Λ0)µ(x, k)eikx dσ. = ∫ ∂Ω ∞∑ l=0 (ik)l l! (Λq − Λ0) ∑ n,m1,m2 ψn,m1,m2e inθkm1k m2 dθ = ∞∑ m1,m2=0 (ik)m1 m1! (ik)m2 m2! tm1,m2 (5.27) Due to the limited domain of convergence of the Taylor expansion, it cannot directly be used to calculate conductivity. Instead, it is necessary to re-express t(k) in terms of a function of the spatial coordinates, c(x). For small perturba- tions in conductivity, ψk(x) ≈ eikx and one may define a linearized conductivity 40 Chapter 5. Approximate Solution (5.28) σ1/2(x) ≈ 1 + c(x), where c' 1. This leads to a useful interpretation of t(k), t(k) ≈ ∫ Ω eikx ∆σ1/2 σ1/2 eikx dx(5.29) ≡ ∫ Ω eikx∆c(x)eikx dx(5.30) = ∫ Ω eikx+ikx4∂∂∆c(x) dx(5.31) By the fundamental theorem of calculus and assuming compact support of c(x) (which is clear for the linearized case, but not so clear for large variations in σ), t(k) = ∫ Ω ( 4∂∂eikx+ikx ) ∆c(x) dx(5.32) =− 4|k|2 ∫ Ω c(x)eikx+ikx dx(5.33) ≡− 4|k|2s(k),(5.34) where s(k) is the ordinary Fourier transform of c(x) evaluated at (−2k1, 2k2). This is closely related to a moment problem for c(x), t(k) =− 4|k|2 ∫ Ω eikx+ikxc(x) d2x =− 4|k|2 ∞∑ m1,m2=0 (ik)m1 m1! (ik)m2 m2! ∫ Ω xm1xm2c(x), d2x (5.35) Recovering c(x) is a linear inverse problem and is solvable with Tikhonov regu- larization. 5.4.1 Low Frequency Approximation The low-frequency approximation for t(k) is t(k) = ∫ ∂Ω eikx (Λq − Λ0) [I + S0(Λq − Λ0)]−1 eikx dσ(5.36) ≡ ∫ ∂Ω eikx+ikx∆c(x)(5.37) 41 Chapter 5. Approximate Solution For any harmonic v and measured boundary potential ψf , (5.38) ∫ ∂Ω v (Λq − Λ0)ψf dσ = ∫ Ω v∆c(x) dx, where u(x) is the harmonic potential satisfying u = f |∂Ω. Each measured volt- age and current pattern may be related to integrals of c(x). Sinusoidal current patterns for radial problems may be used to generate moments of c(x). Then c(x) itself is estimated by solving the inverse problem with Tikhonov regular- ization. Refer to Appendix C for a detailed calculation of the low-frequency approximation of the scattering transform of a disc. 5.5 The ∂-bar Equation 5.5.1 First-order interpretation of µ(x, 0) Once the scattering transform t(k) is found, the next step in reconstruction is to compute µ(x, k) within domain Ω, Equation 4.39 µ(x, k) = 1 + 1 pik ∗ e−k(x)t(k)µ(x, k) 4pik . The convolution in 4.39 may be explicitly written as a double integral, µ(x, k) = 1 + 1 pik ∗ e−k(x)t(k)µ(x, k) 4pik = 1 + 1 (2pi)2 ∫ R2 t(k′) (k − k′) k′ e−x(k ′)µ(x, k′) dk′1dk ′ 2. (5.39) For small perturbations in conductivity, this equation admits a useful inter- pretation in terms of s(k) and c(x). To first order, µ may be set to 1 in the right-hand side, µ(x, k) ≈ 1 + 1 (2pi)2 ∫ R2 t(k′) (k − k′) k′ e−x(k ′) dk′1dk ′ 2 = 1 + 1 (2pi)2 ∫ R2 −4|k|2s(k′) (k − k′) k′ e−x(k ′) dk′1dk ′ 2. (5.40) 42 Chapter 5. Approximate Solution The conductivity is recovered by Equation 4.40, from σ(x) = µ(x, 0)2, i.e. σ1/2 ≈ µ(x, 0) = 1 + 1 (2pi)2 ∫ R2 −4|k|2s(k′) −|k|2 e−x(k ′) dk′1dk ′ 2 =1 + 1 pi2 ∫ R2 s(k′)e−x(k′) dk′1dk ′ 2 =1 + c(x), (5.41) which is consistent with Equation 5.28. 5.5.2 Solution in terms of convolution The following algorithm solves the ∂-bar equation from t(k). Define the Fourier transform pair F [f(x)] = ∫ R2 eikx+ikxf(x) dx ≡ f(k)(5.42) F−1[f(k)] = 1 pi2 ∫ R2 e−ikx−ikxf(k) dk ≡ f(x)(5.43) Then the ∂-bar equation for µ becomes µ(x, k) = 1 + ∂ −1 k ◦ [ t(k) 4pik e−k(x)µ(x, k) ] (5.44) =1 + ∂ −1 k ◦ [−4|k|2s(k) 4pik e−k(x)µ(x, k) ] (5.45) =1 + ∂ −1 k ◦ [−4kF [c(x′)] 4pi e−k(x)µ(x, k) ] (5.46) =1 + ∂ −1 k ◦ [ −iF [∂x′c(x+ x′)]F [µ(x,−x′)] pi ] (5.47) Taking the inverse Fourier transform of both sides, µ(x, y) = δ(y)− 1 pi2 ∫ R2 e−iky−iky∂ −1 k ◦ iF [∂x′c(x+ x′) ∗ µ(x,−x′)] pi dk(5.48) = δ(y)− 1 pi2 ∫ R2 e−iky−iky iF [∂x′c(x+ x′) ∗ µ(x,−x′)] iypi dk(5.49) = δ(y)− 1 piy · [ ∂x′c(x+ x′) ∗ µ(x,−x′) ] (5.50) This last form, Equation 5.50, is implemented numerically. A square grid of size 43 Chapter 5. Approximate Solution 2N is used, with grid spacing h. The convolution is performed as IFFT(FFT(f).∗ FFT(g)) ∗h2, where .∗ is an element-wise product on the grid. The derivative is implemented with central differencing, df/dx1 = [f(x+ h)− f(x− h)]/2h, and the element at the origin of 1/y is set to zero. The conductivity may be recovered from the nonlinear version of 5.41, σ1/2 ≈ µ(x, 0) = 1 + c(x+ x′) ∗ µ(x,−x′) ∣∣∣ x′=0 (5.51) The final computation may be simplified with Parseval’s theorem, however the above form is more useful as a diagnostic. The implementation is sensitive to the interaction between the 1/y term and sharp changes in conductivity. In particular, it appears to be necessary to mollify the singularity in 1/y so that it is smoother than the radius of curvature of discontinuities in conductivity. If not, a Gibb’s-like resonance appears near the edges. Knudsen [22] likely has a superior algorithm. 44 Chapter 5. Approximate Solution Figure 5.3: ∂-bar solution µ(0, x) for a disc object 45 Chapter 5. Approximate Solution 5.5.3 Centre of radially symmetric problems For the centre of radially symmetric problems, t and µ are real and functions of |k| alone [38], and as a result there is an analytical solution. For real µ, Equation 4.39 at the origin becomes (5.52) 1 µ ∂µ ∂k = 1 2µ ( ∂µ ∂k1 + i ∂µ ∂k2 ) = ∂ log |µ| ∂k = t(k) 4pik . Thus the ∂-bar equation for µ is reduced to the first-order equation for for log |µ|. Using the same interpretation for c(x), without the requirement that it be small, but with the requirement that it have compact support, we have at the origin, log |µ(0, 0)| = c(0),(5.53) σ(0) = e2c(0).(5.54) This was tested numerically for a disc shaped object in Figure 5.5.3. Rounding errors limit accuracy for very small conductivities. Also, mollifying the 1/y term in solving the ∂-bar equation tends to reduce µ slightly. 5.6 Conductivity The final step of Nachman’s reconstruction procedure is to calculate the con- ductivity via Equation 4.40, (5.55) σ(x) = µ(x, 0)2. However, an examination of the low-frequency approximation suggests a varia- tion. The scattering transform t(k) is shown in Section to be related to a dis- tribution c(x). For the case of a single disc scatterer, Appendix C analytically computes c(x) to be the disc of magnitude σ−1σ+1 with identical radius. Further- more, by the translation properties of the scattering transform [38], the same is true regardless of the location of the disc. This suggests two approximations. First, one could compute an image di- 46 Chapter 5. Approximate Solution !! !" !# $ # " ! #$ !! #$ !" #$ !# #$ $ #$ # #$ " #$ ! !%$&'()*'+%$&',-.'+/.+012.'3/)+ +%$& ! %$ & 4056./+21'778 9:2+; Figure 5.4: Conductivity at origin computed from ∂-bar equation. 47 Chapter 5. Approximate Solution rectly from c(x) without solving the ∂-bar equation, (5.56) σ(x) = 1 + c(x) 1− c(x) . For small conductivities, σ(x) ≈ 1+2c(x), and for the disc scatterer, −1 < c < 1, so 0 < σ <∞. An improved reconstruction procedure, although with the potential for ad- ditional numerical inaccuracy, follows from the solution to the ∂−bar equation. Section 5.5.3 indicates that for a radially symmetric problem such as a disc, µ(0, 0) = exp(c(0)) and that a disc shaped c(x) becomes a disc shaped µ(x, 0). Thus the approximation could be repeated by substituting logµ(0, 0) for c(x) as before. The new approximation using the ∂-bar equation is (5.57) σ(x) = 1 + logµ(x, 0) 1− logµ(x, 0) . Plots of both approximations are shown in Figures 5.5 and 5.6. The plots are compressed in the vertical conductivity axis to reflect the compression in sensitivity (σ − 1)/(σ + 1). For disc shaped objects in the first plot, the re- construction is shown to be exact (within the limits of the regularization). The annulus, however, exhibits shielding artifacts for very high conductivities. These images are superior to the most recent published work [22]. It is interesting to note that the second approximation with ∂-bar inversion is almost identical to the first, certainly to within the accuracy of the numerical implementation. No significant improvement in image quality was observed as a result of solving the ∂-bar equation for radial conductivities. Perhaps better numerical algorithms would show an improvement, however, it is more likely that the two are equivalent for radial conductivities within the low-frequency approximation of this work. This is interesting because the ∂-bar equation is the slowest part of the reconstruction, as it involves a 2-dimensional calculation for every spatial point in the reconstruction. Instead, the simple approximation of Equation 5.56 provides an accurate one-step linear inverse problem to solve for conductivity, and it is exact for disc shaped objects. 48 Chapter 5. Approximate Solution !! !"#$ " "#$ ! "%%% %%%% "#&% %%%% "#'% %%%% "#(% %%%% "#)% %%%% !%%% &%%% %%%% '%%% %%%% (%%%%%%% %%%%%%%%!"%% &"%% !""" *+,-./012103%45+6178 9:-1.; * + , - . / 01 2 10 3 <=:/0 > " %-!?:5 > " %@AA5+= B1,8:5%@AA5+= (a) !! !"#$ " "#$ ! "%%% %%%% "#&% %%%% "#'% %%%% "#(% %%%% "#)% %%%% !%%% &%%% %%%% '%%% %%%% (%%%%%%% %%%%%%%%!"%% &"%% !""" *+,-./012103%45+6178 9:-1.; * + , - . / 01 2 10 3 <=:/0 > " %-!?:5 > " %@AA5+= B1,8:5%@AA5+= (b) Figure 5.5: Image reconstruction for disc of moderate and high conductivity. 49 Chapter 5. Approximate Solution !! !"#$ " "#$ ! "%%% %%%% "#&% %%%% "#'% %%%% "#(% %%%% "#)% %%%% !%%% &%%% %%%% '%%% %%%% (%%%%%%% %%%%%%%%!"%% &"%% !""" *+,-./012103%45+6178 9:-1.; * + , - . / 01 2 10 3 <=:/0 > " %-!?:5 > " %@AA5+= B1,8:5%@AA5+= (a) !! !"#$ " "#$ ! "%%% %%%% "#&% %%%% "#'% %%%% "#(% %%%% "#)% %%%% !%%% &%%% %%%% '%%% %%%% (%%%%%%% %%%%%%%%!"%% &"%% !""" *+,-./012103%45+6178 9:-1.; * + , - . / 01 2 10 3 <=:/0 > " %-!?:5 > " %@AA5+= B1,8:5%@AA5+= (b) Figure 5.6: Image reconstruction for annulus of moderate and high conductivity. 50 Chapter 6 Conclusion Electrical impedance tomography has evolved considerably from its earliest roots in curve-fitting 1-dimensional electrical sounding data against analytical curves. There has been significant applied research in the geophysical commu- nity for tackling 3-dimensional problems, and new proofs in the mathematical community suggesting strategies for 2 and higher-dimensional problems. New variational methods have been developed, finite-element solutions im- plemented, but even the most modern of these algorithms are essentially a curve fit. Since the problem is extremely ill-posed, information has to be built into the solution method to extract a reasonable image. The approximation presented here reduces 2-dimensional nonlinear inverse conductivity problem to a linear problem, the result of which is curve-fit against an analytic result to construct the image. This approximation is for an idealized circular domain with unit conductivity at the boundary and assumes that arbitrary voltage and current measurements may be made at the perimeter. Even so, the method presented is more general. Much of the approximation would translate well to non-circular or unbounded domains. The final curve-fit assumes constant boundary conductivity, however this too could be generalized. This approximation was motivated from a simple two dimensional test case, a conductive disc embedded at the centre of a circular domain with unit con- ductivity. This problem graphically illustrates the ill-posedness and suggests an image scaling that is later used to map a reconstructed image to conductivity. The approximation evolved from a constructive proof of the uniqueness of the two-dimensional inverse problem by Nachman [33]. Since Nachman’s paper is abbreviated and intended primarily for harmonic analysts, the mechanics of the construction are shown here in detail. This development begins with Calderon’s solution to the linearized inverse conductivity problem. It contains the germ of an idea that motivated further research in this field: by considering combinations of exponential solutions to Laplace’s equation, it is possible to approximate the 51 Chapter 6. Conclusion Fourier transfrom of the conductivity. Beals and Coifman’s inverse scattering theory is then developed, culminating in a formal solution to inverse quantum scattering. Nachman’s construction for conductivity including its relation to boundary data is then shown. Nachman’s result involves exponentially varying functions and a difficult boundary integral. As a result, it is not suitable for direct numerical imple- mentation. An approximation to the construction is developed which allows one to avoid handling exponential boundary voltage patterns, and simplifies the boundary integral. The application of the imaging of a conductive disc results in the scaled image of a conductive disc, regardless of the conductivity of that disc. It is possible to scale the resulting image to reconstruct the conductivity of the original disc. When applied to other shapes, the resulting image is only approximate, but still quite reasonable. Sample reconstructions are illustrated for disc and ring-shaped objects of moderate and high conductivities. The approximation is valid for very high- contrast images, and in principle has no range limit. However, high-contrast objects will tend to shield their interiors. This is intrinsic to the physics of the problem. The resulting images are superior to those presented in the literature for this technique in that they greatly extend the allowable image contrast. The key to achieving this performance is the low-frequency approximation used in the boundary-integral. A significant amount of recent research into this problem has focussed on the latter part of Nachman’s construction, solving the ∂-bar equation. The results in this paper show that, at least for radial conduc- tivities, solving the ∂-bar equation does not materially improve imaging perfor- mance. Rather, improvements to the approximation of the boundary integral are necessary for significant improvement. A side result from developing this approximation is an analytic formula for Fadeev’s Green’s function. Several papers including a chapter of Siltanen’s thesis [38] have been written on calculating this exponentially varying function. The analytic formula presented here relates the Green’s function to the Exponential Integral function, a classical special function whose numerical evaluation has been well studied. 52 Chapter 6. Conclusion 6.1 Summary of the Approximation (a) Measure Dirichlet-Neumann map Λq (6.1) Λq : u|∂Ω → ∂u ∂n ∣∣∣∣ ∂Ω . (b) For each boundary voltage ψf , compute the background potential f sat- isfying (6.2) f = [1 + S0 (Λq − Λ0)]ψf (c) Calculate c(x) satisfying (6.3) ∫ ∂Ω v (Λq − Λ0)ψf dσ = ∫ Ω v∆c(x)u(x) dx. This is a linear inverse problem requiring regularization. (d) Optionally solve for µ(x, y) within domain Ω (6.4) µ(x, y) = δ(y)− 1 piy · [ ∂x′c(x+ x′) ∗ µ(x,−x′) ] . and also µ(x, k = 0) from (6.5) µ(x, k = 0) = 1 + [ c(x+ x′) ∗ µ(x,−x′) ] |′x = 0 or an equivalent expression via Parseval’s theorem. Recompute (6.6) c(x) = logµ(x, k = 0). (e) Calculate conductivity within the domain Ω by scaling based on the ana- lytic expression for a disc, (6.7) σ(x) = 1 + c(x) 1− c(x) . (f) It is suggested that the results be plotted on a scale reflecting the com- pression in the data (σ − 1)/(σ + 1). 53 Chapter 6. Conclusion 6.2 Opportunities for Further Research There are many possibilities for additional research. Extending this approxima- tion to realistic boundary conditions, i.e. using electrode models would be a good step. Afterwards, it should be easy to generalize to unbounded domains. The resolution and sensitivity of the approximations could be studied, however they are sensitive to the method used to solve the linear inverse problem. Anisotropic conductivities and domain errors could also be considered, and perhaps their relationship to each other. The ultimate extension of this approximation would be to apply the technique to Nachman’s 3D inversion algorithm [32], for only in three dimensions would the algorithm become truly useful in applications. 54 Bibliography [1] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Table. Courier Dover Publications, 1965. [2] A. Allers and F. Santosa. Stability and resolution analysis of a linearized problem in electrical impedance tomography. Inverse Problems, 7(4):515– 533, 1991. [3] K. Astala and L. Paivarinta. Calderon’s inverse conductivity problem in the plane. Preprint, 2003. [4] Kendall Atkinson and Weimin Han. Theoretical Numerical Analysis. Texts in Applied Mathematics. Springer, New York, 2001. [5] D.C. Barber and B.H. Brown. Applied potential tomography. Journal of Physics E: Scientific Instruments, 17:723–733, 1984. [6] R. Beals and R. Coifman. Linear spectral problems, non-linear equations and the ∂̄-method. Inverse Problems, 5:87–130, 1989. [7] J.G. Berryman and R.V. Kohn. Variational constraints for electrical- impedance tomography. Physical Review Letters, 65(3):325–328, 1990. [8] B.H. Brown, D.C. Barber, and I.L. Freeston. Tomography, October 21 1986. US Patent 4,617,939. [9] R. M. Brown and G. Uhlmann. Uniqueness in the inverse conductivity problem for nonsmooth conductivities in two dimensions. Communications in Partial Differential Equations, 22:1009–1027, 1997. [10] R.M. Brown, G.A. Uhlmann, R.M. Brown, and G.A. Uhlmann. Uniqueness in the inverse conductivity problem for nonsmooth conductivities in two dimensions. Communications in Partial Differential Equations, 22(5):1009– 1027, 1997. 55 Bibliography [11] A.P. Calderon. On an inverse boundary value problem, Seminar on Numer- ical Analysis and its Applications to Continuum Physics, Soc. Brasileira de Matematica, Rıo de Janeiro, pages 65–73, 1980. [12] S.H. Chan. A study on the direct interpretation of resistivity sounding data measured by Wenner electrode configuration. Geophysical Prospecting, 18(2):215–235, 1970. [13] M. Cheney and D. Isaacson. Distinguishability in impedance imaging. Biomedical Engineering, IEEE Transactions on, 39(8):852–860, 1992. [14] M. Cheney, D. Isaacson, and J.C. Newell. Electrical impedance tomogra- phy. SIAM Rev, 41(1):85–101, 1999. [15] T. Dahlin and B. Zhou. A numerical comparison of 2D resistivity imaging with 10 electrode arrays. Geophysical Prospecting, 52(5):379–398, 2004. [16] John W. Dettman. Applied Complex Variables. Dover, New York, 1965. [17] D.C. Dobson and F. Santosa. An image-enhancement technique for elec- trical impedance tomography. Inverse Problems, 10:317–334, 1994. [18] P. Hua, E.J. Woo, J.G. Webster, W.J. Tompkins, S.G. Inc, and H. Es- tates. Iterative reconstruction methods using regularization and optimal- current patterns in electrical impedance tomography. Medical Imaging, IEEE Transactions on, 10(4):621–628, 1991. [19] D. Isaacson, J.L. Mueller, J.C. Newell, and S. Siltanen. Imaging cardiac activity by the D-bar method for electrical impedance tomography. Physiol. Meas, 27:S43–S50, 2006. [20] K. Knudsen. On the Inverse Conductivity Problem. PhD thesis, Ph. D. thesis, Aalborg University, 2002. [21] K. Knudsen. A new direct method for reconstructing isotropic conductivi- ties in the plane. Physiological Measurement, 24(2):391–401, 2003. [22] K. Knudsen, M. Lassas, J.L. Mueller, and S. Siltanen. D-bar method for electrical impedance tomography with discontinuous conductivities. Tech- nical report, Dept. of Mathematical Sciences, Aalborg University, 2006. 56 Bibliography [23] K. Knudsen, J. Mueller, and S. Siltanen. Numerical solution method for the dbar-equation in the plane. Journal of Computational Physics, 198(2):500– 517, 2004. [24] O. Koefoed. A generalized cagniard graph for the interpretation of geoelec- trical sounding data. Geophysical Prospecting, 8(3):459–469, 1960. [25] O. Koefoed. A semi-direct method of interpreting resistivity observations. Geophysical Prospecting, 13(2):259–282, 1965. [26] O. Koefoed. The direct interpretation of resistivity observations made with a wenner electrode configuration. Geophysical Prospecting, 14(1):71–79, 1966. [27] R. Kohn and A. McKenney. Numerical implementation of a variational method for electrical impedance tomography. Inverse Problems, 6(3):389– 414, 1990. [28] R.E. Langer. An inverse problem in differential equations. Bull. Amer. Math. Soc, 39:814–820, 1933. [29] M.H. Loke and R.D. Barker. Rapid least-squares inversion of apparent re- sistivity pseudosections by a quasi-Newton method 1. Geophysical Prospect- ing, 44(1):131–152, 1996. [30] J.L. Mueller and D. Isaacson. Regularization of the computed scattering transform for the D-bar method for electrical impedance tomography. Pro- ceedings of SPIE, 5562:121, 2004. [31] J.L. Mueller, S. Siltanen, and D. Isaacson. A direct reconstruction algo- rithm for electrical impedance tomography. Medical Imaging, IEEE Trans- actions on, 21(6):555–559, 2002. [32] A. Nachman. Reconstructions from boundary measurements. Ann. Math, 128(2):531–576, 1988. [33] A. Nachman. Global uniqueness for a two-dimensional inverse boundary value problem. Annals of Math, 143:71–96, 1996. [34] K. Paulson, W. Breckon, and M. Pidcock. Electrode modelling in electrical impedance tomography. SIAM J. Appl. Math, 52(4):1012–1022, 1992. 57 Bibliography [35] N. Polydorides and W.R.B. Lionheart. A Matlab toolkit for three- dimensional electrical impedance tomography: a contribution to the Electri- cal Impedance and Diffuse Optical Reconstruction Software project. Mea- surement Science and Technology, 13(12):1871–1883, 2002. [36] Walter Rudin. Real and Complex Analysis. McGraw-Hill, New York, 1987. [37] S. Siltanen, J. Mueller, and D. Isaacson. An implementation of the recon- struction algorithm of A. Nachman for the 2D inverse conductivity problem. Inverse Problems, 16(3):681–699, 2000. [38] Samuli Siltanen. Electrical Impedance Tomography and Fadeev Green’s Functions. PhD thesis, Helsinki University of Technology, 1999. [39] S. Stefanesco and C. Schlumberger. Sur la distribution électrique potentielle autour d’une prise de terre ponctuelle dans un terrain à couches horizontales homogènes et isotropes. Journ. de Phys. et du Radium, VII:132–140, 1930. [40] A.N. Tikhonov. On the uniqueness of the solution of the problem of electric prospecting. Dokl. Akad. Nauk SSSR, 69:797–800, 1949. [41] E.J. Woo, P. Hua, J.G. Webster, andW.J. Tompkins. A robust image recon- struction algorithm and its parallelimplementation in electrical impedance tomography. Medical Imaging, IEEE Transactions on, 12(2):137–146, 1993. 58 Part II Appendices 59 Appendix A Schrödinger equation and Laplace’s equation A.1 Theorem. Let σ(x) be twice differentiable. The Laplace equation with spa- tially varying conductivity, (A.1) ∇ · (σ∇v) = 0, may be transformed into an equivalent Schrödinger equation, (−∆+ q)u = 0, where u = σ1/2v and q = ∆σ1/2/σ1/2. Proof. This is actually easier to see in reverse. Substitute u = σ1/2v in (4.21). ∆u− qu = ∆(σ1/2v)− ∇σ 1/2 σ1/2 σ1/2v = ∇ · ( σ1/2∇v + vσ −1/2 2 ∇σ ) − v∇ · ( σ−1/2 2 ∇σ ) = σ1/2δv + σ−1/2 2 ∇σ ·∇v + σ−1/2 2 ∇v ·∇σ − vσ −3/2 4 ∇σ ·∇σ + vσ −1/2 2 ∆σ + v σ−3/2 4 ∇ ·∇ σ − vσ −1/2 2 ∆σ = σ1/2∆v + σ−1/2∇σ ·∇v = σ−1/2(σ∆v +∇σ ·∇v) = σ−1/2(∇ · (σ∇v)). 60 Appendix B Potential for One and Two Discs Single Disc An analytical expansion for the potential of a disc will be developed in terms of a general solution for piecewise-constant, radially symmetric regions. As a starting point, it is useful to calculate all possible solutions to a potential equation with the test object embedded in an infinite 2-dimensional medium of unit conductivity. This is most easily done by finding all solutions asymptotic to xn as x→∞. A useful starting point is the Kelvin transform. B.1 Definition. The Kelvin Transform maps the interior of the unit disc to its exterior, and vice versa, T (x) = x|x2| , x = (x1, x2) The transform may be generalized to invert about a disc centred at xa with radius ra, Ta(x) = raT ( x− xa ra ) = r2a x− xa |x− xa|2 + xa In the notation of complex variables, with {x ≡ (x1, x2)}→{ z ≡ x1 + ix2}, Ta(z) = r2a 1 z − za + za. If f(x) is harmonic on the interior of a disc, the Kelvin transform f(T (x) is harmonic on the exterior of the disc. This elementary result may be found in Atkinson [4]. For a conductive disc in a planet of uniform conductivity, the potential will be harmonic on both the interior and exterior of the disc. A perfect conductor will 61 Appendix B. Potential for One and Two Discs have constant potential on its interior (du/dn = 0 on the interior), a perfect insulator will have no normal electric field at its surface (du/dn = 0 on the exterior), and a constant conductivity will have a discontinuity in du/dn based on the contrast with the surrounding medium. Without the disc, the potential solutions for homogeneous 2-dimensional medium are entire functions. With a conducting disc, one may express the potential as a sum of an entire function and a decaying portion due to the disc. A useful function to describe this is given as follows B.2 Definition. Let v be a potential harmonic inside a disc Da, where Da is centred at xa and has radius ra. The scattering primitive is defined by w(x, v) = { (v(x)− v(xa)), |x− xa| < ra (v(y)− v(xa)), |x− xa| > ra, , where y = Ta(x). The following properties make the function useful for solving potential prob- lems with discs. B.3 Theorem. The scattering primitive, w, is harmonic in the interior and exterior of the disc, continuous across the disc boundary, and satisfies ∂w ∂n ∣∣∣∣ ∂Ω+ = − ∂w ∂n ∣∣∣∣ ∂Ω− , Proof. ∂w(y) ∂n(y) = ∂w (y(x)) ∂n (y(x)) = ∂w(x) ∂n(x) ∂n(y) ∂n(x) = −∂w(x) ∂n(x) by the inversion property of the Kelvin transform. Given an entire function v then v±w is a solution to a potential problem in the interior and exterior of a unit disc. The properties of w are such that v−w has constant potential interior to the disc, and v+w has zero normal derivative, i.e. these are solutions to the conductivity problem for a perfect conductor and perfect insulator. This may be further generalized to constant conductivities. B.4 Theorem (Single disc scattering, low-frequency approximation). Let a disc of conductivity σa and radius ra be placed at position xa in plane of conductivity sigma0. The true potential satisfies the laplace equation ∆ · (σ∆v) = 0. If the 62 Appendix B. Potential for One and Two Discs background potential is vb, then the scattered potential is fw(x, vb) where f is a material parameter, f = 1− σaσ0 1 + σaσ0 . Proof. At the disc boundary, the Laplace equation reduces to v|∂Ω+ = v|∂Ω− , and σ0 ∂v ∂n ∣∣∣∣ ∂Ω+ = σa ∂v ∂n ∣∣∣∣ ∂Ω− , where ∂Ω+ and ∂Ω− denote the interior and exterior of the disc. The background potential and its normal derivative are continuous across the disc boundary. The function w is continous, but its first derivative switches sign. Inside the disc, the first derivative of w coincides with vb. The linear combination of potentials vb+ cw asymptotically approaches vb, and is continuous across the interface. We also have, ∂(vb + fw) ∂n ∣∣∣∣ ∂Ω− = (1 + f) ∂vb ∂n ∣∣∣∣ ∂Ω− and ∂(vb + fw) ∂n ∣∣∣∣ ∂Ω+ = (1− f) ∂vb ∂n ∣∣∣∣ ∂Ω− . Therefore, 1 + f 1− f = σa σ0 , or f = 1− σaσ0 1 + σaσ0 . Note that a disc of infinite conductivity will “cancel” its background potential (f = −1), while a disc of zero conductivity will “reflect” its background potential (f = 1). Figure B.1 shows scattering from a metallic disc. Two Disc Scattering The scattered potential from two discs will be developed in terms of a series of single-disc images. Each disc will independently respond to the background po- tential, then to the scattered potential from the other disc, then to the scattered- scattered potential, and so forth. The result is a summable series for a scattered potential consistent with the boundary conditions and background potential. 63 Appendix B. Potential for One and Two Discs (a) Background potential (b) Scattered potential (c) True potential Figure B.1: A harmonic potential is scattered from a metallic disc. The ap- plied, scattered, and true potentials are shown in B.1(a), B.1(b), and B.1(c) respectively. The metallic disc is an equipotential, and the scattered potential decays. 64 Appendix B. Potential for One and Two Discs One may picture a train of waves hitting two rocks in a pond. The ripples are the sum of the initial train striking the rocks and from the reflected waves slosh- ing between them. Representative background, scattered, and total potentials are shown for two metallic discs in Figure B.3. B.5 Theorem. Consider two disjoint discs Dk, k = {1, 2} with radii rk and centres xk. The Kelvin transform of D2 in D1 is a disc D ⊂ D1. Proof. Since the Kelvin transorm is the conjugate of a linear fractional transformation, L(z) = a+ bz c+ dz , the proof is nearly identical to the corresponding proof in complex variables, see Dettman [16]. It will provide a useful estimate. Consider the family of lines and circles in the plane. α(x2 + y2) + βx+ γy + δ = 0. or in the notation of complex variables, α |z|2 + β ( z + z 2 ) + γ ( z − z 2i ) + δ = 0. Under the transformation w = T (z) = 1/z, the curves become α 1 |w|2 + β 2 ( w + w |w|2 ) + γ 2i ( w − w |w|2 ) + δ = 0. δ |w|2 + β ( w + w 2 ) + γ ( w − w 2i ) + α = 0. δ(u2 + v2) + βu+ γv + α = 0, If the disc a disc is located outside the unit disc, α and δ will be non-zero, so the transformation maps the exterior disc to a disc on the interior. The general Kelvin transform is just a shifted and scaled verision of w = 1/z. Figure B.2 shows two disjoint discs, their Kelvin transformed images, and higher order images. B.6 Theorem. If the separation of D2 and D1 is d, the radius of D is r = r1r22 (r2 + d)(r2 + 2r1 + d) . 65 Appendix B. Potential for One and Two Discs Figure B.2: Reflected images of two discs. and r r2 ≤ 1 2 . Proof. Construct a unit vector between the centres of the two discs. ξ̂ = x2 − x1 |x2 − x1| . The nearest and farthest points ofD1 are mapped to two points ofD interior to D2. Using a geometrical approach, x1 + r1ξ̂ → x2 − r2 r2 r2 + d ξ̂ x1 − r1ξ̂ → x2 − r2 r2 r2 + d+ 2r1 ξ̂. So the radius of D is r = 1 2 [ r2 r2 r2 + d − r2 r2 r2 + d+ 2r1 ] = r1r22 (r2 + d)(r2 + d+ 2r1) . B.7 Theorem. Consider a disc D3 ⊆ D1. Let D3 be reflected intoD2 to produce D4 and then back into D1 to produce D. Then the radius of D satisfies r r3 ≤ ( r1 r1 + d )2( r2 r2 + d )2 , 66 Appendix B. Potential for One and Two Discs where d is the separation between the two discs. Proof. This estimate arises by applying B.6 twice, taking r3 to be small, and noting that the separation between both D3 and D2, and D4 and D1 is greater than d. The blue discs of Figure B.2 illustrate how successive reflections shrink in size. B.8 Theorem. Let Wv = w1 (·, w2(·, v)). If v is harmonic in D1, then sup R2 |Wnv| ≤ mαn, n = 1, 2, . . . for some constants m > 0 and 0 < α < 1. Proof. Let m = 1 2r1 sup D1 |∇v| = 1 2r1 sup ∂Ω1 ∣∣∣∣ ∂v∂n ∣∣∣∣ . then for any two points x and x′, sup D1 |v(x)− v(x′)| ≤ m 2r1 |x− x′| ≤ m. By successive applications of Theorem B.7, Wn(x, v) maps the R2 into the sequence of range v(Bn) where Bn ⊆ D1, and the radii are constrained by rBn r1 ≤ [( r1 r1 + d )2( r2 r2 + d )2]n . Setting α = ( r1 r1 + d )2( r2 r2 + d )2 < α < 1, we obtain sup D1 |Wn(x, v)| = supx,x′∈Bn |v(x)− v(x′)| =≤ 2rn m2r1 |v(x)− v(x ′)| ≤ mαn. B.9 Theorem (Two disc scattering). Consider two discs, Dk k ∈ {1, 2}, placed in a plane of conductivity σ0. The scattered potential may be developed via a series of interactions between the two discs. 67 Appendix B. Potential for One and Two Discs Proof. Let k′ denote the index of the other disk, i.e. if k = 2, k′ = 1. Let vk denote the potential scattered by Dk. By Theorem B.4, (B.1) vk = fkwk(·, vb + vk′). The total potential is (B.2) v = vb + v1 + v2. We expand vk in a sum of multiple reflections vik such that v0k = 0 vi+1k = fkwk(·, vb + vik′). If vik → v̂k uniformly, then v̂k will satisfy B.1, i.e. v̂k = vk. We may reforormulate the vk expression as v0k = 0 v1k = fkwk(·, vb) vi+2k = fkfk′wk (·, vb + wk′(·, vik)) ≡ fkfk′wk ◦ wk′(·, vb + vik). This may be reformulated as by separating out each reflection. Define rik so that (B.3) vik = i∑ j=0 rik. Then r0k = 0 r1k = fkwk(·, vb) ri+2k = fkfk′wk (·, wk′(·, vb + rik)) = fkfk′W (·, vb + rik). 68 Appendix B. Potential for One and Two Discs (a) Background potential (b) Scattered potential (c) True potential Figure B.3: A harmonic potential is scattered from a two discs. The applied, scattered, and true potentials are shown in B.3(a), B.3(b), and B.3(c) respec- tively. or r2nk = (fkfk′) nWn(·, vb), n = 1, 2, . . . r2n+1k = (fkfk′) nfkW n(·, wk(·, vb)), n = 1, 2, . . . By Theorem B.8, rk → 0 uniformly in R2 with geometric convergence. Hence vik → vk uniformly. Furthermore, by Harnack’s Theorem [36], vk is harmonic both within Dk and outside of Dk. 69 Appendix C Single disc scattering transform, low-frequency approximation C.1 Theorem. Let a disc of radius r1 and conductivity σ1 be placed at the the origin in a medium of conductivity 1. The scattering transform of the disk is, t(k) = 4pifr1 |k|J1(2r1 |k|), where f = ( 1−σ1 1+σ1 ) . Proof. If the background potential is eikx, the corresponding scattering potential is (C.1) w = eik r21 x̄ . t(k) = ∫ ∂Ω eik̄x̄ (Λq − Λ0)µ(x, k)eikx ds = ∫ ∂Ω eik̄x̄ (Λq − Λ0) ( eikx + feikr 2 1/x̄ ) ds = ∫ ∂Ω eik̄x̄ (Λq − Λ0) ( eikx + feikr 2 1/x̄ ) ds = ∫ ∂Ω ∞∑ m=0 ( (ikx)m m! ×−2|n|f ∞∑ n=1 (ikr21x)n n! ) ds =− 4pif ∞∑ n=1 (−1)nkr12n n!(n− 1)! 70 Appendix C. Single disc scattering transform, low-frequency approximation Using the definition of the Bessel function J1, (C.2) −z 2 J1(z) = ∞∑ n=1 (−1)n ( 12z)2n−1 n!(n− 1)! , from which we get (C.3) t(k) = 4pif |k|r1J1(2|kr1|). According to the interpretation in Equation 5.29, this may be related to a function c(x) and its transform s(k), (C.4) s(k) = − t(k) 4|k|2 = − pifr1J1(2|k|r1) |k| . The fourier transform of a disc is related to J1, (C.5) ∫ ra 0 ∫ 2pi 0 ei2k1r cos θ dr dθ = piraJ1(2k1ra) k1 , which when combined with the radial symmetry of the disc at the origin shows that c(x) is a disc with radius ra and magnitude −f , i.e. (C.6) c(x) = σ1 − 1 σ1 + 1 |x| ≤ r1 0 |x| > r1 Note that for σ1 = 1 + δσ1, δσ1 ' 1, (C.7) 1 + c(x) ≈ { 1 + δσ1/2 ≈ σ1/21 |x| ≤ r1 1 |x| > r1 Furthermore, Siltanen shows in [38] that for two conductivities related by translation, σ2(x) = σ1(x− x′), the scattering transforms shift as would a linear fourier transform, t2(k) = ek(x′)t1(k). Therefore the low-frequency scattering transform for a disc, regardless of position, is a scaled and trans- lated version of the linearized scattering transform for a disc. 71
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- An approximation method for electrical impedance tomography
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
An approximation method for electrical impedance tomography Pereira, Paulo J. S. 2008
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | An approximation method for electrical impedance tomography |
Creator |
Pereira, Paulo J. S. |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | Electrical impedance tomography is an imaging method with applications to geophysics and medical imaging. A new approximation is presented based on Nachman's 2-dimensional construction for closed domains. It improves upon existing approximations by extending the range of application from resolving 2 times the surface conductivity to imaging perfect conductors and insulators. With perfect knowledge of boundary data, this approximation exactly resolves a single conductive disc embedded in a homogenous domain. The problem, however, is ill-posed, and imaging performance degrades quickly as the distance from the boundary increases. The key to the approximation lies in (a) approximating Fadeev's Green's function (b) pre-processing measured voltages based on a boundary-integral equation (c) solving a linearized inverse problem (d) solving a d-bar equation, and (e) scaling the resulting image based on analytical results for a disc. In the development of the approximation, a new formula for Fadeev's Green's function is presented in terms of the Exponential Integral function. Also, new comparisons are made between reconstructions with and without solving the d-bar equation, showing that the added computational expense of solving the d-bar equation is not justified for radial problems. There is no discernible improvement in image quality. As a result, the approximation converts the inverse conductivity problem into a novel one-step linear problem with pre-conditioning of boundary data and scaling of the resulting image. Several extensions to this work are possible. The approximation is implemented for a circular domain with unit conductivity near the boundary, and extensions to other domains, bounded and unbounded should be possible, with non-constant conductivity near the boundary requiring further approximation. Ultimately, further research is required to ascertain whether it is possible to extend these techniques to imaging problems in three dimensions. |
Extent | 1551939 bytes |
Subject |
Electrical impedance Electrical resistance Tomography |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-08-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066563 |
URI | http://hdl.handle.net/2429/1536 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2008-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2008_fall_pereira_paulo.pdf [ 1.48MB ]
- Metadata
- JSON: 24-1.0066563.json
- JSON-LD: 24-1.0066563-ld.json
- RDF/XML (Pretty): 24-1.0066563-rdf.xml
- RDF/JSON: 24-1.0066563-rdf.json
- Turtle: 24-1.0066563-turtle.txt
- N-Triples: 24-1.0066563-rdf-ntriples.txt
- Original Record: 24-1.0066563-source.json
- Full Text
- 24-1.0066563-fulltext.txt
- Citation
- 24-1.0066563.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066563/manifest