Accurate D-bar Reconstructions of Conductivity Images Based on a Method of Moment with Sinc Basis.

Planar D-bar integral equation is one of the inverse scattering solution methods for complex problems including inverse conductivity considered in applications such as Electrical impedance tomography (EIT). Recently two different methodologies are considered for the numerical solution of D-bar integrals equation, namely product integrals and multigrid. The first one involves high computational burden and the other one suffers from low convergence rate (CR). In this paper, a novel high speed moment method based using the sinc basis is introduced to solve the two-dimensional D-bar integral equation. In this method, all functions within D-bar integral equation are first expanded using the sinc basis functions. Then, the orthogonal properties of their products dissolve the integral operator of the D-bar equation and results a discrete convolution equation. That is, the new moment method leads to the equation solution without direct computation of the D-bar integral. The resulted discrete convolution equation maybe adapted to a suitable structure to be solved using fast Fourier transform. This allows us to reduce the order of computational complexity to as low as O (N (2)log N). Simulation results on solving D-bar equations arising in EIT problem show that the proposed method is accurate with an ultra-linear CR.

INTRODUCTION

Numerical solution to scattering equation is the key in solving a variety of inverse problems in engineering and science including electromagnetic inverse scattering,[1] quantum inverse scattering,[2345] and medical imaging.[6] There are many medical imaging applications including cross-well tomography,[7] ground penetrating radar surface imaging,[89] ultrasound tomography,[10] optical microscopy,[1112] photon imaging,[13] induction imaging,[14] positron emission tomography (PET)[15] and electrical impedance tomography (EIT).[6161718192021]

There are several numerical methods for solving non-linear inverse scattering problem including direct and iterative approaches. D-bar is one of the significant methods based on direct methodology introduced by Beals and Coifman[22] that recently received a lot of attentions due to circumventing highly complex iterative approaches. D-bar is based on formulating several non-linear inverse scattering integrals into some linear D-bar equations.[2324]

Planar D-bar equation in partial differential form can be represented by:[22]

Where the solution μ: R2 → C to the D-bar equation is asymptotic to 1. That is, when k is large, μ(k) ≈ 1. Here T (k) is a complex function. The -operator in (1.1) is defined by:

The variable k may be defined in the real domain (k ∈ R2) or complex domain (k ∈ C); in other words in Eq. 1.2, the notations k = (k1, k2) and k = k1 + ik2 are interchangeable. Convolving both sides of Eq. 1.1 with Green's function for -operator[22] and using asymptoticity of μ the integral equation of D.bar becomes:

where k = k1 + ik2 ∈ C,k ≠ 0

A couple of methods are introduced for the numerical solution of the D-bar equation including product integrals (PI) and multigrid (MG). The PI method is introduced by Siltanen et al.,[25] first, for implementing an EIT algorithm based on the constructive proof of Nachman[24] for two-dimensional inverse conductivity problem. In this approach, a uniform grid is used to discretize the Eq. 1.3. Then, it computes the singular convolution integral in D-bar equation via separating it into smooth and singular parts. The singular part of the integral is computed analytically and the smooth part may be computed using an interpolator polynomial. Then the discrete form of the equation may be solved by means of an iterative solver such as generalized minimal residual solver (GMRES).[26] Note that, the PI based methods requires O (N6) arithmetic computations to converge to the approximate solution of the D-bar equation on N-point grids, which is huge even for advanced ultra-fast computers.

In order to circumvent the huge computational complexity of PI methods, more recently, another class of algorithms is developed by Knudsen et al.[27] based on multigrid method. Note that, MG methodology was first introduced by Fedorenko as claimed[28] in 1962. Derivations of this methodology have shown good efficiency in solving partial differential equations (PDEs).[2930] In the year 2000, over three decades later, Vainikko used MG methods for solving Lippmann-Schwinger integral equations.[31] In 2004, Knudsen et al. adapted Vainikko's MG method for solving the D-bar equation as mentioned earlier. The first step in MG method is to obtain the periodic version of the D-bar integral equation. Next, the periodic equation is discretized using a special kind of grid mapping in which each function in (1.3) is expanded in terms of Fourier basis on the points of a uniform grid. This type of mapping allows neglecting the singular computations of the convolution kernel; however, it imposes some error in approximating the solution of (1.3). Although fast implementations of MG approach gain a remarkable speed and decrease the computational burden from O (N6) to O (N4log N) via the use of fast Fourier transforms (FFT) on N-point grids, these methods suffer from low convergence rate (CR), especially near discontinuities. For example, the CR of the adaptation of MG in[27] is of O (h) where h denotes the spacing parameter of the grid, which is considered very low.

Disadvantages of the MG method motivated us for a new computational method. We found out that sinc-based methods could guarantee ultra-linear CR s in solving integral equations.[323334353637] As a result, we decided to use the method of moments with sinc basis to solve the planar D-bar integral Eq. 1.3.

Therefore, in this paper a method of moments with sinc basis functions speed moment method (SMOM) is used to numerically solve the D-bar Eq. 1.3. This is based on sinc interpolation function; and it is able to achieve an ultra linear CR.[38] The sinc based methods are already used for solving specific ordinary differential equations, PDE's and integral equations in diverse applications in science and engineering.[353637383940414243] Sinc methods show considerable efficiency in terms of accuracy and computational complexity especially in tackling singular computations.[38]

In the proposed method, the nonpareil features of sinc functions are considered to solve the planar D-bar equation, efficiently. We show that the orthogonality of sinc basis functions may lead to discretize Eq. 1.3. Thus, the closed form solution of (1.3) is not required and so, it can help us to circumvent the direct computations of two-dimensional D-bar integral equations. In the following section, the details of the proposed algorithm is introduced and discussed.

METHOD OF MOMENTS WITH SINC BASIS FUNCTIONS

Disceretizing the D-bar Equation

In this section, we derive a solution for μ(k) from Eq. 1.3 based on the method of moments. We start by expanding the product of . Assume the product is frequency band.limited and can be expanded in terms of a set of sinc basis function as

where[373844]

Similarly, the Green's function can be expanded as

Here, h is the constant step size between grid points in each coordinates kq, q = 1, 2. Inserting expansion (2.1) in Eq. 1.3 gives us:

Next, substituting (2.3) in (2.4) gives

Remarkably, the orthogonality of sinc basis[3738] helps us to get

Substituting Eq. 2.6 in Eq. 2.5 leads to

Finally, Eq. 2.7 can be evaluated at discrete grid points k = (mh, nh) as

The above equation is a discrete two-dimensional convolution equation. In the following section, we consider the computational complexity of solving Eq. 2.8 assuming that T (k) = 0 outside an open disk B (0, R). We show that this assumption helps to reduce the computational complexity of the solution.

Reducing the Computational Complexity

Let us denote the domain of two-dimensional integral in Eq. 1.3 by Ω and support of T (k) by supp (T). Unlike the continuous function g (k), the samples of T (k) are usually available via some computations on experimental data in a finite number of points k, where k ∈ supp (T). For example, in the D-bar integral equation for EIT, T (k) may be obtained via finite samples in k space using “discrete voltage to the current map.”[1945]

In addition, due to the commutative property of the convolution integral,[4647] we can convolve T (k) μ (k) over the kernel, instead. That is, we may re-write the Eq. 1.3 as:

where k = k1 + ik2 ∈ C,k ≠ 0

Note that the solution μ (k) for k ∈ supp (T) suffice. That is, the integral Eq. 2.9 may be solved without the need for the values of μ (k) at points outside the problem bound, i.e., k ∉ supp (T). As a result, all we need to compute the convolution (2.9) is to supply the discrete values of the product T (k) μ (k) at points k ∈ supp (T).

In addition, the integral equation may be solved easier in the polar coordinate system using the notation k = k1 + ik2 instead of the Cartesian notation k = (k1, k2) in Eq. 2.9 as well as in Eq. 1.3. Now, based on our own crucial observations, if sup (T) is assumed to be embedded in an open disc B (0, R) (i.e., supp (T) ⊂ B (0, R)), then the convolution integral of (2.9) does not require the values of the kernel for |k'| ≥ 2R. herefore, we may only need to find convolution values inside the open disk of radius 2R as it is shown in Figure 1. More specifically, as it is illustrated in Figure 1b, to compute the convolution integral in the right hand side of the Eq. 2.1, there is no need to use the values of the kernel function outside of the open disc B (0, 2R). As a result, we characterize the domain of integration as Ω = [–2R, 2R] × [–2R, 2R] and rewrite Eq. 1.3 as

Figure 1

Bounds of two-dimensional convolution in D-bar equation and required grid structure. (a) Domain of integration and support of T (k). (b) T (k) *g (k) in the support of T (k). (c) Required grid

where k = k1 + ik2 ∈ C,k ≠ 0

and Eq. 2.9 as

where k = k1 + ik2 ∈ C,k ≠ 0

That is, the computation grid is finite and as illustrated in Figure 1c may be embedded in the open disk B (0, 2R). The consequence is that the two-dimensional discrete convolution operations in the right hand side of (2.8) may be implemented using discrete Fourier transform (DFT).[4647]

Since the Eq. 2.8 is a discrete two-dimensional convolution equation, the FFT and its inverse (FFT-1) may be used to speed up the DFT computations. Thus, the Eq. 2.8 may be implemented as

Traditionally, the FFT algorithms are more effective when they are dealing with matrices that contain 2m data points in each dimension.[46] Thus in the SMOM procedure, the number of points in the grid, N, is chosen to be powers of two to accelerate FFT computations. On the other hand, extending the number of data points in each dimension of the matrix of the two-dimensional functions, T (k) and g (k), via zero padding may result in a more approximate computation of their DFTs.[46] Therefore, as shown in Figure 2, before doing FFT computations of Eq. 2.12, the related matrices are extended via padding at least N data points in each dimension. The number of double precision mathematical operations required to compute each two-dimensional FFT-based convolution computation is O (N2log N).[46] In addition, FFT (g) is only computed once and reused in all iterations. Moreover, to solve Eq. 2.12 using an iterative solver like GMRES, one should set the upper bound of the number of iterations as a parameter. For example, in GMRES, this parameter should be set as an integer in the range [1…10]. Therefore, only O (N2log N) double-precision mathematical operations are required in each of the iterations of the Eq. 2.12 of the SMOM method, using FFT and FFT inverse. That is, solving Eq. 2.12 for the value of μ (k) in each of N2 discrete points k = k1 + ik2 where 1 ≤ k1, k2≤ N, using an iterative solver with a fixed upper bound on the number of iterations requires O (N2log N) arithmetical operations. It is clear that solving Eq. 2.12 for all N2 discrete points requires O (N2) iteration of the solver. Therefore, FFT is called O (N2) times. As a result, entirely O (N4log N) arithmetical operations are required for N2-point SMOM solution to D-bar equation.

Figure 2

Illustration of zero padding before fast Fourier transforms based implementation of the speed moment method. (a) The proper padding for a function (b) The correct result convolution

It is worth noting that, since FFT algorithms subdivide the computational data into small segments and combine them hierarchically, the effect of round-off error in their computations is considerably minimized.[46]

EIT USING SMOM

In this section, the proposed SMOM method is employed to solve the inverse conductivity problem. The aim of this problem is to reconstruct the bounded conductivity distribution γ of a body Ω via the knowledge of injected current patterns and resulted voltages on the electrodes surrounding the body.[21] An in-depth review for the mathematical formulation of two-dimensional inverse conductivity EIT problem may be found in.[171821] We refer the readers to these sources for more details.

Briefly, let us suppose a bounded domain Ω ⊂ R2. Let also apply a voltage potential on the boundary of Ω, through electrodes of an EIT measurement system. The consequent of this system is a potential distribution u ∈ H1 induced in Ω, which can be obtained as the unique solution of the conductivity equation

Another effect is the normal current distribution which can be measured through the same set of the electrodes. Using the distribution of the potential and normal current on the boundary of Ω, the Dirichlet-to-Neumann map Λ may be defined as

The forward and inverse problems of EIT are related to this mapping. The forward EIT problem is to determine Λ from the knowledge of γ in the mapping γ → Λ; while the inverse problem is to determine γ from the knowledge of Λ in the inverse mapping Λ → γ. The inverse problem was first stated by Calderon.[48]

The significance of the forward EIT problem comes into account when a numerical test model is necessary for the solution evaluation of an inverse problem. In such a case, the Eq. 3.1 is solved as the model solution, with a pre-assumed distribution γ to find the normal current distribution g and normal potential distribution f on the boundary of Ω. Then, any intended inverse conductivity solution may be compared with this model.

For the inverse solution, Brown and Uhlman[49] proved that it is possible to solve such problems with a constructive solution and offered an algorithm. In the following, we review their constructive proof for inverse conductivity problem in C1 conductivities. In this type of conductivity distribution, the induced potential is modeled with a continuous function.

In Brown and Uhlmann algorithm, two D-bar equations arises on the inverse conductivity problem defined on a body. There are some efforts to solve these two equations including Knudsen et al. and Siltanen et al. methods.[2725] In this section, we employ and implement a new algorithm based on SMOM approach for solving the Brown and Uhlmann described version of the two-dimensional D-bar equations.

In the following, we first review Brown and Uhlmann constructive description of the two dimension D-bar problem. Then a numerical model solution implying the cross-section of human chest during expiration, which is suggested in,[27] is formed to evaluate recent D-bar solution methods as well as our SMOM method and compare them together.

Brown and Uhlmann EIT Method

The method of Brown and Uhlmann solves the inverse conductivity problem in two constructive steps:[49]

In the first step, the problem is transformed to inverse scattering problem to arrive at an intermediate function S. In the second step, the conductivity distribution γ is computed from S. Note that the ill-posedness of EIT problem comes from the first step. For the purpose of comparison, we adapted some notations and formulations from.[27]

Let u be a solution to the conductivity equation in the plane. In the following, we identify (z1, z2) ∈ Ω with the complex number z = z1 + iz2. Defining the potential

the is the solution to the system.[49]

A helpful assumption in this method is taking γ = 1 outside the boundary of Ω With this assumption, the potential q can be extended to whole R2 by q = 0, in R2\Ω The scattering transform S is related to potential q via[49]

where . The solutions to the first D-bar equation[49]

are related to the functions m1 (z, k) eikz and m2 (z, k) eikz via[49]

These exponentially growing solutions, with asymptotic conditions[49]

are introduced by Faddeev.[50] Using Eq. 3.8, m1 (z, k) may be computed as

In order to compute the conductivity γ the second D-bar equation is to be solved that may be represented as:[49]

Finally, the unique solution m + to this equation is used to compute[27]

Next, a particular conductivity model of chest cross-section during expiration is defined and our SMOM solution is offered for the related D-bar equations.

Test Example: Numerically Simulated Phantom Chest

Our test example is a numerically simulated cross-section of the phantom chest in a unit circle that is illustrated in the Figure 3a. This numerical phantom is used recently in some EIT researches.[2751] To have a simulated model of the chest conductivity during expiration, the conductivity values of the heart, lungs and background tissue are assumed to be equal to the values of the second column of Table 1. These conductivities are then scaled by dividing by the value of the background conductivity. The scaled conductivity values are enlisted in the third column of Table 1.

Figure 3

The reconstructed test conductivity using speed moment method: (a) Original conductivity γtrue (b) Conductivity γSMOM

Table 1

The conductivity values of the simulated chest phantom

The lungs are modeled by ellipses and the heart is modeled by a circle. For this conductivity the corresponding potential q is computed using numerical differentiation on a 512 × 512 grid defined on z variable. To compute S (k), the forward solution, on the points of the grid which is defined on k variable, first, the Eq. 3.7 is solved on the aforementioned grid defined on z variable. Note that 512 × 512 fine grid guarantees a high accuracy solution to Eq. 3.7. The integral form of Eq. 3.7 is[27]

Next, we compute the function m1(z, k) using Eq. 3.10 and then, evaluate integral (3.6) to compute S (k) using trapezoid rule.

In order to solve the inverse problem and to compute the conductivity in the second step of the reconstruction algorithm m+(z, k) is computed for each z in a new uniform grid in z-plane. Note that, using these two different grids on z variable mitigates the well-known inverse crime problem.[5253] Having computed the scattering transform S (k), the integral equation

is solved for each z in the uniform grid and the conductivities γrec are computed via evaluation of the solution of (3.14) to (3.12) at k = 0.

In the next section, first the implementation of SMOM algorithm as well as PI and MG methods are described. Next, using the subsequent numerical results derived in the application of all three methods on the aforementioned test problem, some significant parameters of efficiency are computed for the proposed SMOM method as well as other two methods and then compared.

SMOM IMPLEMENTATION AND EVALUATION

In order to implement and evaluate our SMOM method for solving D-bar equations in the inverse conductivity problem and compare it with best of its competitor methods, namely PI and MG methods, three Matlab codes were developed by authors one for each. Note that to implement a Matlab code for the PI numerical method, the analytical computation of 16 singular integrals in the form:[25]

Is required. Theses integrals have integrable singularity when the variable s lies in one corner of the mesh element that is defined by [ui, ui + 1] × [νi, νi + 1]. A simple approach that may facilitate to analytically compute these singular integrals is to divide the domain of integration in to two parts as it is done in Figure 4. Computing integrals (4.1) over each In, n = 1, 2 and adding the results gives us the value of (4.1) over the mesh element. Table 2 enlists the results of analytical computation of the integrals (4.1) where the parameters (a, b) = (0, 0) and the evaluation variable s lies in four corners of the mesh element. The value of these singular integrals are analytically computed for other 12 cases and enlisted in Tables 3–5.

Figure 4

Dividing the domain of integration to compute and evaluate the singular integrals (4.1) on the corners of mesh elements in the product integrals method

Table 2

The result of analytical computation of singular integrals

Table 3

The result of analytical computation of singular integrals

Table 5

The result of analytical computation of singular integrals

Implementation

In order to implement and later compare, the iteration parameters should be set in a common range for all three methods namely PI, MG, and our SMOM method. For this purpose, same-size grids are used. In addition, the upper bound of the iterations and acceptable tolerance of the iterative solver was set to 10 and 10−7, respectively, for all three cases.

In order to reconstruct the test conductivity, we follow the steps of the Brown and Uhlmann that explained in the previous section. That is, after constructing the required grid structure, all three methods solve the first D-bar Eq. 3.13. Then, the scattering transform is computed according to the Eq. 3.6 using trapz function of Matlab. Note that Eq. 3.6 can be written as[54]

In this equation, Λ, Λl are the Dirichlet-to-Neumann maps of respectively the inhomogeneous chest phantom with conductivity distribution γ and a homogenous phantom with constant conductivity γ = 1. Employing Eq. 3.6 to compute the scattering transform allows us to circumvent the explicit computation of the Dirichlet-to-Neumann mapsΛ, Λl.

Having computed this intermediate function, the second D-bar equation which is represented by Eq. 3.14 is solved for points k ∈ [-25,25]×[-25,25]. Finally, evaluating the solutions of the second D-bar equation to Eq. 3.12 at k = 0 results the image of the conductivity values in positions of the uniform grids defined on z-space. The result of the reconstruction of the test conductivity via SMOM using 64 × 64 gird in z-plane is illustrated in Figure 3b.

Evaluation

In order to evaluate the proposed method and weigh it against its competitor methods, the most important and frequently-used parameters of efficiency[55] namely, accuracy, dynamic range of reconstructed conductivity (DRRC) and computational complexity versus accuracy are examined.

Accuracy

Due to the outperformance of MG method compared with PI methodology, we only discuss and compare the accuracy of SMOM with that of MG method. Similar to MG method, we follow the Knudsen et al. methodology,[27] for SMOM accuracy. That is, we use the supermom norm as a comparative measure of accuracy for the solution of the D-bar equations in the aforementioned EIT test problem as was used in the MG methods. In addition, since the true value of m+(z, k) is not available for all k, and also to avoid the quantization errors that emerges in the computation of S (k), the solution matrix m+(.,1) is considered in both methods. That is, first, the matrix m+(.,1) is computed as the solution to Eq. 3.14 at k = 1 with a 512 × 512 grid in k-plane. Next, the solution matrix m+(.,1) is used as the true solution against the other approximate solutions that are resulted with coarser grids in z-plane. Let us denote this true solution as mtrue(.,1) and denote the approximate solution with a grid of size Ni, i = 1,…,5 as mi(.,1). The supermom norm of solution error may be defined as:

In addition, the CR can be defined as:

Thus, the Eq. 3.14 is solved at k = 1 using the SMOM with different levels of disceretization in k-plane that are enlisted in the second column of the Table 6. For each solution, the value of Ei is computed as it is reported in the third column of the Table 6. According to the Eq. 4.4, the values of error norm are used to compute the respective CRs that are shown in the fourth column of Table 6. The results show that the SMOM follows an ultra-linear CR.

Table 6

Analysis of the convergence rate of the approximate solution m (.,1) computed by the SMOM method to the true solution m (.,1)

Table 7 shows the analysis of the CR of MG method using supermom norm reported by Knudsen et al.[27] Comparing the Table 7 results with that of the SMOM in Table 6, we observe that the new proposed SMOM method possesses a higher CR than MG.

Table 7

Analysis of the convergence rate of the approximate solution m (.,1) computed by the one-grid method to the true solution m (.,1)[27]

Another important parameter to evaluate and rank an EIT numerical solution method is the analysis of the DRRC of resulted images. Such analysis is offered in many EIT publications including[16] and.[45] In fact, DRRC can reflect the efficiency of the numerical solution method used for solving the D-bar equations. In the following subsection, we shall use this parameter and analyze our SMOM as well as PI and MG numerical methods, and compare them together.

DRRC

Let define the DRRC as:[1645]

This parameter demonstrates the potency of the reconstruction algorithm in separating the different levels of the reconstructed conductivity. The diagram in the Figure 5 compares the resulted values of DRRC parameter in different reconstructions of all three methods at different levels of disceretization in k -plane. Computing the DRRC of SMOM, PI and MG methods, it can be observed that the DRRC using the new SMOM method is closer to 100% compared to each of the other two methods. For example, while the DRRC of the PI and MG with 128 × 128 grids in k-plane are about 74% and 92.2% respectively, the DDRC of SMOM is about 95.56%.

Figure 5

Dynamic range of the reconstructed conductivities via applying the product integrals, multigrid and sinc based moment method

Computational complexity

For comparing the computational complexity of our method with that of PI and MG, we measured the required computing time (in seconds) for solving Eq. 3.13 at a fixed point k = 1 using a dual core, 2.53 GHz processor speed Sony laptop. For this purpose, the following scenario is arranged:

Set a common environment for all three methods. That is, equal numbers of grid points are set. In addition, the acceptable tolerance and the upper bound of the iterations of the solver are set to 10−6 and 10, respectively

Under common settings of part (1), D-bar equation defined in (3.13) is solved at point k = 1 with different levels of disceretization.

All three methods' performances are reported in Table 8. The results indicate that the computational complexity of the new SMOM is way too smaller than that of the PI method. In addition, SMOM takes much less computing time than MG method, as well.

Table 8

The time required for computing the approximate solution of equation (3.13) in the fixed point k = 1 using different sizes of grids

CONCLUSION

Applying the -operator on non-linear equations, which arise in quantum inverse scattering, PET and EIT -result into a convolution form D-bar integral equation. High computational complexity of the PI methods and low precision of the MG methods demonstrate the inefficiency of both of them in solving the D-bar equation.

In this work, sinc functions are utilized as basic functions to introduce a new sinc basis Method of Moment (SMOM) for solving the planar D-bar equation. The brilliant feature of this novel method is that, upon the orthogonality of sinc basis functions, the computation of highly computational-burden convolution integral in the D-bar equation is circumvented. In this method, the discrete two-dimensional convolution equation of D-bar easily lends itself to a fast computational implementation via the FFT.

The performance of our method is compared with the best available methods namely PI and MG. In order to numerically compare these methods, a scenario was adapted for our SMOM algorithm as well as the other two methods, similar to settings adapted by other researchers. That is, we adapt an inverse conductivity problem defined on a numerically simulated chest phantom.

Numerical results confirm the effectiveness and the superiority of our SMOM method for solving planar D-bar equations in terms of both the computational complexity and accuracy over MG and PI.

BIOGRAPHY

Mahdi Abbasi is an assistant professor of Bu-Ali Sina University, Hamedan, Iran. He received the B.Sc. degree in computer engineering and M.Sc. in computer Architecture from Sharif University of Technology, Iran, in 2000 and 2005 respectively. He received Ph.D. in computer Architecture, from University of Isfahan, in 2012. His areas of research interest include medical imaging systems and algorithms, image processing, network processor systems and algorithms.

E-mail: abbasi@basu.ac.ir

Table 4

The result of analytical computation of singular integrals