Next Article in Journal
Evolutionary Algorithm with Geometrical Heuristics for Solving the Close Enough Traveling Salesman Problem: Application to the Trajectory Planning of an Unmanned Aerial Vehicle
Previous Article in Journal
A Pixel-Wise k-Immediate Neighbour-Based Image Analysis Approach for Identifying Rock Pores and Fractures from Grayscale Image Samples
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The D-Bar Algorithm Fusing Electrical Impedance Tomography with A Priori Radar Data: A Hands-On Analysis

1
Helmholtz Institute for Biomedical Engineering, RWTH Aachen University, 52074 Aachen, Germany
2
Department of Physics, Goethe University Frankfurt, 60438 Frankfurt, Germany
*
Author to whom correspondence should be addressed.
Algorithms 2023, 16(1), 43; https://doi.org/10.3390/a16010043
Submission received: 30 November 2022 / Revised: 29 December 2022 / Accepted: 5 January 2023 / Published: 9 January 2023

Abstract

:
Electrical impedance tomography (EIT) is an imaging modality that can estimate a visualization of the conductivity distribution inside the human body. However, the spatial resolution of EIT is limited because measurements are sensitive to noise. We investigate a technique to incorporate a priori information into the EIT reconstructions of the D-Bar algorithm. Our paper aims to help engineers understand the behavior of the D-Bar algorithm and its implementation. The a priori information is provided by a radar setup and a one-dimensional reconstruction of the radar data. The EIT reconstruction is carried out with a D-Bar algorithm. An intermediate step in the D-Bar algorithm is the scattering transform. The a priori information is added in this exact step to increase the spatial resolution of the reconstruction. As the D-Bar algorithm is widely used in the mathematical community and thus far has limited usage in the engineering domain, we also aim to explain the implementation of the algorithm and give an intuitive understanding where possible. Different parameters of the reconstruction algorithm are analyzed systematically with the help of the GREIT figures of merit. Even a limited one-dimensional a priori information can increase the reconstruction quality considerably. Artifacts from noisy EIT measurements are reduced. However, the selection of the amount of a priori information and the estimation of its value can worsen the reconstruction results again.

1. Introduction

Electrical impedance tomography (EIT) is a modality that allows estimating the conductivity distribution inside the human body. Compared to other established imaging techniques, such as magnetic resonance imaging or computed tomography, it has some distinct advantages: firstly, an EIT system can be constructed at a low cost; secondly, it does not depend on potentially harmful X-rays; thirdly, it can be applied long-term at the bedside of the patient [1]. Medical applications of EIT may include the monitoring of the lung, including observing regional lung ventilation and perfusion [2], estimation of the size and volume of the bladder [3], monitoring of lung recruitment and lung collapse [4], tracking of 3D brain activity [5], and usage in breast cancer imaging [6]. The EIT can also be used outside the medical context, for example, for monitoring semiconductor manufacturing [7] or monitoring damage in concrete [8].
The EIT belongs to the class of inverse problems. The solution to such problems requires unique algorithms. These methodologies can be grouped into three major groups. The first set of methods are the variational regularization algorithms. These methods minimize a cost function that is made up of two parts. The first part minimizes the given measurements to a conductivity distribution, while the second part is a regularization term. Two examples of variational regularization algorithms are Thikonov regularization [9] and the total variational regularization [10]. The variational methods are mostly used in the field, as they deliver good results in rather quick time. Commercial devices such as the Draeger PulmoVista® 500 use this kind of algorithm. The second group is the statistical inversion methods for EIT. The measurements and conductivities are modeled as random variables. Here, an a posteriori distribution can be estimated, which, in turn, can give estimates of the conductivity distribution [11]. The advantage of these types of algorithm is that not only is a reconstruction generated, but so is a whole statistical distribution of the reconstruction. Some publications also generate better results with such methods [12]. However, a major drawback is the computational cost. The third type of reconstruction algorithms are the direct inversion methods. They analyze the nonlinear inverse problem and try to develop an algorithm based on this analysis. One such algorithm is the D-Bar algorithm [13].
The D-Bar method is a direct, noniterative algorithm used to compute the conductivity distribution inside a domain of interest. The insights and tools for this algorithm date back to Calderón [14], who analyzed the inverse conductivity problem. Subsequently, Nachmann showed that the problem is uniquely solvable in the two-dimensional (2D) domain under certain restrictions [15]. The first implementation of this algorithm was shown by Siltanen et al. and, thus, opened the door for real-world usage of the algorithm [16]. Nowadays, the D-Bar algorithm is known to have many benefits, such as robustness towards electrode displacement, measurement errors, and errors in the injection current [17].
Significant disadvantages of EIT are its relatively low spatial resolution and sensitivity to noise. Regularization strategies are used to dampen the noise dependence. The regularization is a trade-off between a further loss of spatial resolution and resistance to noise. Due to this fact, there are several strategies to incorporate information from other modalities. This is often carried out by biasing the solution process with the a priori information. Sun et al. modified the sensitivity matrix by using X-ray data and a collection of impedance measurements to improve the reconstruction results [18]. Adler and Lionheart proposed nudging the solution of the inverse problem towards the desired solution by careful selection of the prior [19]. This approach could also craft a custom regularization matrix [18]. Soleimani proposed a level set method based on ultrasound as a priori data [20]. Alsaker and Mueller proposed an a priori information strategy based on the D-Bar algorithm from CT data [21]. The stimulation patterns in the D-Bar algorithm used are the trigonometric patterns; these are selected because of the orthogonality to each stimulation. Although this is motivated on a theoretical level, it is confirmed by the results of Adler et al. on the work of optimal injection patterns [22,23].
In this paper, we want to explain the D-Bar algorithm through a hands-on engineering approach and analyze the behavior of one-dimensional a priori data in a 2D reconstruction domain. We take a simplified version of their proposed algorithm and analyze its behavior based on additional radar data, which adds further information along the x-axis. We explain the implementation of a D-Bar algorithm and sweep across some parameters to gain a more intuitive understanding of the algorithm. Since the radar data in our case only affect a priori information along the x-axis, the behavior of the D-Bar algorithm exhibits interesting properties. In the case of breast cancer, radar data provide more detail and are less noisy; however, compared to EIT, the ease of usage is more complicated and damping of radar waves in fat tissue is an issue. As radar data still provide valuable additional information, we used the setup from Patri et al. for our a priori data [24]. We note that Flores-Tapia and Pistorius have already embedded radar information into EIT reconstructions, but without using the D-Bar algorithm [25].

2. Fundamentals

2.1. Basic EIT Equations

The EIT handles the problem of reconstructing the conductivities inside a domain of interest, for example, the human body, through boundary voltage measurements. Accordingly, N electrodes are placed equidistantly across a slice of a body. A current pattern is injected through these electrodes, and the resulting voltages are measured. There are different current injection and voltage measurement patterns, which have distinct advantages and disadvantages. For more information on this, see [26]. The voltage measurement patterns also have different advantages and disadvantages. The problem of EIT can be modeled as follows:
( σ ( z ) u ( z ) ) = 0 , z Ω
where u represents the electrical potential, Ω is the 2D domain (the slice through the body) of interest, σ is the location-dependent conductivity of the domain, and z denotes the position expressed as a complex number z = x + j · y (with j = 1 ). The position is denoted as a complex number z due to mathematical considerations with the usage of complex geometrical optics solutions. Only real-valued conductivities are considered for the sake of simplicity. The system from Equation (1) is excited by an electric current to reconstruct an estimate of σ inside Ω . Mathematically speaking, this is the Neumann boundary condition that constrains the partial differential equation of Equation (1). This can be expressed in the notation of Equation (1) as
σ ( z ) u n ( z ) = J ( z ) , z Ω
where n is the outward normal vector of the boundary of Ω , J is the current density, and Ω denotes the boundary of Ω . In theory, it is also possible to apply known voltage distributions on the boundary of Ω and then measure the resulting currents. The known voltage distribution is called the Dirichlet boundary condition and can be expressed in the notation of Equation (1) as
u ( z ) = f ( z ) , z Ω
where f ( z ) is an arbitrary function in the sense that, depending on the injected current, it can lead to different forms. It is possible to map the Dirichlet boundary condition to the Neumann boundary condition. This map is called the “Dirichlet-to-Neumann map” (DN-map). It can more formally be stated as
Λ σ : u σ u n
Astala and Päivärintas showed that knowledge of the DN-map dictates the conductivity inside Ω uniquely [27]. Our take on the explanation of its fundamentals can be found in the Appendix A.1 or, for a more detailed look, in Chapters 14 and 15 of Mueller’s and Siltanen’s book [28].

2.2. D-Bar Algorithm

The D-Bar algorithm can be summarized into four steps:
  • Compute the DN map Λ σ from the voltage measurements.
  • Compute the scattering transform within | k | R 1 .
  • Solve the D-Bar equation, Equation (A4).
  • Recover the impedance values from Equation (A6).

3. Numerical Implementation

3.1. Numerical Implementation of the D-Bar Algorithm

We now describe the numerical implementation of the D-Bar algorithm used in this work. All simulations were carried out using EIDORS [19]. We provide the original functions in Matlab® used for this paper as Supplementary Materials along with this paper.
Throughout this paper, N denotes the number of electrodes, while M denotes the number of linear independent stimulation current patterns. In the setting of EIT, M = N 1 holds true as the maximum number of linear independent current patterns. Note that this is true for the chosen trigonometric current pattern. Different choices of patterns may result in a reduced number of linear independent patterns. An example of a stimulation pattern is given in Figure 1 with a c o s i n e stimulation around a circular domain. The current through each electrode is given by c o s ( 8 · Θ n ) (n denotes the n t h electrode), which is similar to the case m = N / 2 in Equation (5) (m denotes the m t h stimulation pattern). The standard circular domain was used for better visualization.
Step 1: Compute the DN map Λ σ from the voltage measurements and known stimulation currents. Λ σ is an operator that can map a given set of boundary voltages to its resulting boundary currents. Λ σ is approximated through a matrix L σ with the dimension N 1 × N 1 [28]. Thus, to map a discrete voltage measurement to its corresponding current stimulation, the vector of voltage measurements on each electrode can be multiplied by L σ to obtain the corresponding current distribution. We first compute the Neumman-to-Dirichlet (ND) matrix R σ and then invert it to obtain L σ , because the direct computation of the DN matrix L σ involves derivatives. Let J be the matrix of stimulated current distributions, whose columns contain one stimulation pattern each with the appropriate current on every electrode. J is defined as
J n m = { I cos ( m Θ n ) , m = 1 , , N / 2 1 I cos ( π n ) , m = N / 2 I sin ( ( m N / 2 ) Θ n ) , m = N / 2 + 1 , , N 1 ,
where m is the column index (one specific stimulation), n is the row index (the electrode position), I is the current amplitude, and Θ n is the position angle of the n t h electrode from the model center. Since the maximum number of linear independent measurements is M = N 1 and the number of electrodes is N, the matrix has a shape of N × N 1 = N × M .
The matrix of voltages measured U describes the voltage on every electrode for each stimulation pattern and has the same shape as J . Each column of U contains the voltages measured on the electrodes. The voltage on the n t h electrode for the m t h stimulation pattern is u n m . Regarding the m t h stimulation pattern, the voltages on the electrodes must fulfill n = 1 N u n m = 0 . U can be assembled with [ U ] n , m = u n m .
We also need a geometrical correction term (see Figure 2), which scales each entry of R σ with the arc length differential, since our domain is not circular. This technique was introduced by Herrera et al. [29]. The arc length differential for each electrode position n and current injection m is calculated by
s n m = r 2 ( Θ n ) + r ( Θ n ) r ( Θ n 1 ) Δ Θ n , n 1 2 · Δ Θ n , n 1 r 2 ( Θ n ) + r ( Θ n ) Θ n 2 · d Θ ,
where r ( Θ n ) is the radius (distance from model center to electrode) of the model at electrode angle Θ n . The correction matrix S can be assembled by [ S ] n , m = s n m . As can be seen in Equation (6), s n m is independent of m. This is because the geometry does not change with each current stimulation. The results are then formed to a vector, which is repeated M times to form the N × M matrix U .
R σ can be calculated with the help of J , U , and S through
R σ = 1 A ( J S ) T U .
A denotes the electrode area, while ⊙ is the elementwise multiplication and T represents the transpose.
The last step is to invert R σ to finally obtain the DN matrix L σ = R σ 1 . The calculation steps for the DN matrix can be efficiently performed in Matlab® through matrix multiplication.
Step 2: Calculation of the scattering transform.
The scattering transform can be viewed as a complex Fourier transformation. Instead of Equation (A5), we use an approximation of the scattering transform t e x p ( k ) = Ω e i k ¯ z ¯ ( Λ σ Λ 1 ) e i k z d s , which can be calculated according to Dodd and Mueller [30]: The functions e i k ¯ z ¯ | Ω and e i k z | Ω need to be expanded in their orthonormalized current pattern basis. The coefficients used are c ( k ) = [ c 1 ( k ) , , c N ( k ) ] T for e i k z and d ( k ) = [ d 1 ( k ) , , d N ( k ) ] T for e i k ¯ z ¯ . The coefficients are calculated with
c ( k ) = J e i k z a n d d ( k ) = J e i k ¯ z ¯ .
The exponential approximation of the scattering transform is then expressed by
t e x p ( k ) i = 1 N 1 j = 1 N 1 c j ( k ) d i ( k ) ( L σ ( i , j ) L 1 ( i , j ) )
= d T ( k ) L σ L 1 c ( k ) .
The scattering transform represents the data in a frequency-like complex-valued space called the k-space. Unlike the Fourier transformation, which is linear, in the sense that it represents a function through sinusoidal waves, the scattering transform represents the input via nonlinear waves. t e x p ( k ) is only evaluated for values of | k | < R , where R is called the truncation radius. This is due to the fact that high values of | k | lead to numerical instabilities. An example of the k-space for a specific conductivity enclosure and | k | < 3.5 can be observed in Figure 3.
For the implementation of this computational step, a grid in the k-space has to be implemented. This becomes important, as choosing the grid wisely helps in the computation of the next step. In this work, we chose the strategy of Mueller and Siltanen, which is presented in Chapter 15 of their book [28], for the selection of the grid. This strategy was used in this work. In the next step of the D-Bar algorithm, a singularity arises in the origin in the k-space. For the computations inside the k-space, a matrix K is formed. It has the shape P × P = 2 p × 2 p . An example is visualized in Figure 4. This has the advantage that the singularity at the origin is avoided. The values inside K are computed depending on the chosen truncation radius R. This is carried out with the help of the Matlab® function meshgrid(). This function obtains a vector k h e l p = [ R , R + 2 R M 1 , , R 2 R M 1 , R ] as an input argument. Finally, we can construct K = K 1 + j · K 2 with K 1 , K 2 = m e s h g r i d ( k h e l p ) . Note that K is complex-valued; this is for computational purposes, as it makes the following steps easier. A simplified example of K can be seen in the Appendix A.2.
Step 3: Solution of the D-Bar Equation (A4).
As mentioned in the Appendix A, the D-Bar algorithm solves the conductivity for one position. Thus, the steps outlined in the following have to be repeated for every position desired inside the domain. We use the method proposed by Mueller and Siltanen, which is based on a fast Fourier transform approach [28]. Equation (A4) can be regarded as the convolution of the differential form of Equation (12) and the Green’s function G ¯ for the D-Bar operator, leading to
μ ( k , z ) = 1 + G ¯ ( T ( k , z ) μ ( k , z ) ¯ ) . w i t h
T ( k , z ) = t ( k ) 4 π k ¯ e z ( k )
G ¯ = 1 π k
Note that the x–y-position inside the domain is expressed as a complex number, with z = x + j · y . Instead of regarding Equation (11) for the full plane, it is transformed to a periodic equation. To achieve this, a grid Q = [ 2 R 3 ϵ , 2 R + 3 ϵ ] 2 with ϵ > 0 is defined. This ensures a sufficient grid size for the calculations to perform the FFT trick by Mueller and Siltanen [28]. In our case, ϵ = R / 6 was used. Q is very similar to the earlier constructed K , except that it is larger.
Q = 2 R 3 ϵ + j · ( 2 R + 3 ϵ ) 2 R + 3 ϵ + j · ( 3 R + 3 ϵ ) 2 R 3 ϵ j · ( 3 R + 3 ϵ ) 2 R + 3 ϵ j · ( 3 R + 3 ϵ ) .
Since the fast Fourier transformation is used later, Q should be of the shape 2 p × 2 p , where p can be freely chosen. We used p = 6 for this work, as it resembles a good compromise between accuracy and computation speed. In addition, note that the grid—similar to K —is symmetricand, thus, the origin 0 + j 0 is not in the grid. This simplifies calculations, otherwise, the origin is a singularity in the calculations (see Equations (12) and (13)).
On this grid, the Green’s function (13) is evaluated elementwise. We denote the resulting matrix as G ¯ to make the function represented by the matrix periodic. The matrix is weighted such that the boundary of the matrix is zero, while a circle inside the matrix has the original values. The weight function η ( k ) is defined as
η ( k ) = { 1 | k | < 2 R + ϵ ( 2 R + 3 ϵ | k | ) / ( 2 ϵ ) 2 R + ϵ | k | . 0 2 R + 3 ϵ < | k |
The periodic version of G ¯ is denoted as G ˜ ¯ and calculated through G ˜ ¯ ( k ) = G ¯ ( k ) η ( k ) , where ⊙ denotes the elementwise multiplication (“.*”—operation in Matlab®).
Now that the periodic version of Equation (13) G ˜ ¯ ( k ) is calculated, the other known part of Equation (11) T ( k , z ) needs to be evaluated. Note that K is a subset of the larger grid Q . For every grid point of Q , T ( k , z ) = t ( k ) 4 π k ¯ e z ( k ) is evaluated. Note that we still write T ( k , z ) rather than using Q. This is because T ( k , z ) is only different from 0 inside the truncation radius. Thus, we use K for the computation and just pad the border with 0 to make the dimensions compatible with the rest of the operations on Q . This can be performed through the Matlab® function padarray. As explained earlier, this step (step 3) has to be repeated for every position z, where we want to know the conductivity. Thus, there is a small gain to make by computing the first part of T ( k , z ) ( T f i r s t ( k ) = t ( k ) 4 π k ¯ ) beforehand and saving it to memory.
Now, the FFT trick from Mueller and Siltanen can be applied on Equation (11) using the fast Fourier transform, through
μ = I + h 2 I F F T ( F F T ( G ˜ ¯ ) F F T ( T μ ¯ ) )
I =   I μ P ( T μ ¯ ) w i t h
P ( T μ ¯ ) = h 2 I F F T ( F F T ( G ˜ ¯ ) F F T ( T μ ¯ ) ) ,
where h is the spacing of the grid K , I is the identity matrix, and ⊙ denotes the elementwise multiplication similar to that in Equation (7).
This equation can be solved with the help of a matrix-free equation solver such as GMRES, as the function cannot be brought to the simple A x = b form. This is due to the fact that Equation (17) has a complex conjugate of μ in it. Thus, the equation is only linear in its real part. When using GMRES, the real and imaginary parts must be detached from one another. This can be achieved by transforming μ from its matrix form to its vector form and separating the imaginary part and the real part. The complex matrix μ with dimensions C H × H will become a vector of dimension R ( 2 · H 2 ) , and the vector representation is expressed as μ . The first M 2 entries are the real components of μ , while the last M 2 are its imaginary components. The mapping from the vector form to the matrix form is denoted as v e c m a t . This procedure is also described in [28].
A normal GMRES call in Matlab® to solve A x = b would result in gmres(A, b). This call must be modified. Instead of a matrix A, a function f is given to gmres(). The function has two tasks; firstly, it converts μ from its vector to the matrix representation μ ( v e c m a t ); secondly, it calculates the first step inside the GMRES loop, which is called the “Krylov vector”. The calculation is given by
q n e w = q o l d v e c m a t ( P ( T q o l d ¯ ) ) ,
where q is the “Krylov vector” from the GMRES method. For more details about GMRES, see [31].
Step 4: Calculation of the conductivity.
This is performed through σ ( z ) = μ 2 ( k 0 , z ) . This can be understood when observing Equation (A2). Ψ ( k , z ) = μ ( k , z ) holds true for k = 0 , then the equation becomes ( Δ + σ 1 / 2 Δ σ 1 / 2 ) μ ( k , z ) = 0 . Since σ ( z ) = μ ( z , 0 ) , the equation can be reformulated to
( Δ + Δ σ ( z ) σ ( z ) ) σ ( z ) = 0 Δ σ ( z ) + Δ σ ( z ) = 0 .
Thus, a suitable solution to (A2) is σ ( z ) = μ ( z , 0 ) . Squaring both sides of the equation leads to the calculation of σ ( z ) .
The process of solving the D-Bar equation needs to be repeated for every position z inside the domain which will be reconstructed.

3.2. EIT Simulation Setup

The motivation behind the fusion is to be able to visualize tumors in the female breast. Radar, as additional information for EIT, can increase the resolution of EIT. Relying solely on radar, however, might not be ideal as fat dampens the radar signals. The general setup without any tumor inclusions for this paper is equal to the model shown in Figure 5. The electrodes are placed in the cyan-colored areas. To mimic the female breast, as in Prati et al. [24], we also added a 2 mm skin layer. The inside of the breast was modeled as fat tissue. The conductivity values for both tissues were extracted from the Gabriel database [32]. The conductivity enclosure is a tumor, with the values from the measurements of Surowiec et al. [33].

3.3. Radar Simulation Setup

For the simulation of radar, the setup from Prati et al. was used [24]. The antenna patches are squared and have an edge length of l p a t c h = 3 mm. The antennas are etched on the substrate RO4350B with a dielectric constant of ϵ r = 3.55 @ 10 GHz. The antenna is excited by a coaxial probe, whose position was chosen such that the input impedance matched at 50 Ω . Radar antennas were placed on the top and bottom of the domain Ω . The bottom antennas are used as transmitters, and the antennas on the top are used as receivers. A pulse excitation was used for the stimulation of the antennas. One excitation cycle was composed of a sawtooth signal with an up-chirp time of 433 μ s , a down-chirp time of 100 μ s , and a waiting time in between the chirps of 100 μ s . The simulation was carried out as follows.
The described radar antennas were only used on the bottom and top of the domain. While a fixed distance between the radar antennas is required for the above-described radar setup, in reality, it is more practical to produce a circuit board with embedded antennas already on it. Since breasts have different sizes, the placing of radar antennas on the sides of the domain requires strenuous effort, to still ensure the correct placement of the antennas. EIT electrodes, on the other hand, can be placed manually on the sides, as they can be designed to stick on the skin.
As in the EIT setup, the skin was modeled with a thickness of 2 mm. For the skin, the complex dielectric properties of the Gabriel database were used [32]. The breast fat and tumor tissues were modeled from the values of Martellosio et al. of the subgroup with high adipose tissue [34]. The radar simulations were performed in CST Microwave Studio 2019.

3.4. Radar Reconstruction

At first, a reference simulation with no enclosed tumor was performed as described above. Each transmission antenna sends a pulse signal through the medium without inclusion. The signals were captured by the receivers on the top. The reference signals are denoted as s a , b 0 = [ s 1 0 , s 2 0 , , s T 0 ] , where a indicates which transmitter was used for sending the signal, b denotes the transmitter receiving the signal, and T denotes the length of the signal.
After that, the simulation was carried out with the inclusions. In our case, we investigated one inclusion. The inclusion signals are denoted as s a , b = [ s 1 , s 2 , , s T ] , using the same notation as in the first case.
The reconstruction process was carried out as follows. At first, the root-mean-square deviation (RMSD) was computed for each transmitter–receiver pair ( a , b ) by
R M S D a , b = 1 T t = 1 T ( s a , b 0 ( t ) s a , b ( t ) ) 2 .
The vector R M S D a is defined as R M S D a = [ R M S D a , 1 , , R M S D a , B ] , where B is equal to the number of receivers. The next step for the reconstruction is to compute the mean R M S D ¯ a = 1 B b = 1 B R M S D a , b of each R M S D a for each transmitter, where B is the number of receivers. Only values above the mean are used for the final reconstruction, resulting in
R M S D a , b f i n a l = { R M S D a , b R M S D a , b > R M S D ¯ a 0 o t h e r w i s e .
The cumulative R M S D a , c u m for every receiver position is calculated from R M S D a , b f i n a l with R M S D a , c u m = b = 1 B R M S D a , b f i n a l . The position of the enclosure can be estimated from the resulting curve, as it shows local minima. This can be seen in Figure 6. An envelope of the RMSD curve is constructed to obtain a better estimation of the tumor positions. This is performed by gathering the local maxima and interpolating linearly between them. The position estimate of the tumor is given by subtraction of the envelope with the RMSD curve. This estimation is used as a priori information in the D-Bar algorithm.

3.5. Fusion of EIT and Radar

The spatial resolution of the radar is higher than that of EIT due to the low wavelength used by radar. However, only information about the position of the object along the x-axis is given in the introduced reconstruction of radar data, while the EIT is able to capture the spatial details in both axes. Thus, a careful weighing of EIT and radar information is necessary. A fusion can be performed by usage of the k-space and the scattering transform t ( k ) . Regarding | k | R 1 , we calculate the scattering transform for the EIT t R 1 E I T ( k ) using Equation (9). Regarding R 1 < | k | R 2 , we calculate t R 2 r a d a r ( k ) . For that, the conductivity inside the domain needs to be estimated. This was performed as follows. The peak of the radar reconstruction (see red line in Figure 6) gives the most likely position of the tumor; moving farther off the peak, the tumor position is less likely (the curve is decreasing). Along the x-axis of the FEM phantom, we added the difference conductivity of tumor and fat to the domain weighted by the R M S D a , b f i n a l -curve. Along the y-axis, the values are just replicated, or, in other words, along the y-axis, values are padded such that vertically the values of the conductivity do not change. From a strict mathematical point of view, we have a two-dimensional prior; however, as the information is fixed along one axis, for all practical purposes, information is only carried in one dimension.
After that, the calculations follow the procedure described in Section 3. This mean that we use Equation (10) to calculate t e x p ( k ) . Thus, the calculation of the prior is performed the same way as for the EIT data. The new scattering transformation with EIT and radar information can be stated as
t f u s i o n ( k ) = { t R 1 E I T ( k ) , | k | R 1 t R 2 r a d a r ( k ) , R 1 < | k | R 2 0 , R 2 < | k | .
t f u s i o n ( k ) combines the advantages of both the EIT and radar. The concept is also visualized in Figure 7. Note that this is equivalent to Alsaker et al.’s Equation 3.6 in [21]. However, the actual calculation of t R 2 r a d a r ( k ) is different from their calculations. In addition, we have a one-dimensional prior instead of a two-dimensional one, compared to their setup. Later, we analyze the behavior of under- and overestimation of the conductivity, by scaling the added conductivity with a factor. The reconstruction with the help of t f u s i o n follows the same procedure as the pure EIT reconstruction (step 3 of the D-Bar algorithm).

3.6. Evaluation of the Reconstructions

The fusion of D-Bar and EIT data with the D-Bar algorithm introduces a variety of hyperparameters that all affect the quality of the reconstruction. The standard D-Bar method has the truncation radius R 1 as a hyperparameter. In our approach, a second truncation radius is introduced, which is similar to [35]. As mentioned in Section 3.5, the a priori conductivity σ p r i o r also has an effect on the reconstruction. Finally, the effect of noise on the reconstruction is of great interest, as the radar data might be able to make the algorithm more robust to noise.
We use the GREIT figures of merit for the evaluation of the reconstructions [36]. These performance metrics are well established for the evaluation of different reconstruction algorithms in EIT. The data have to be preprocessed to use the GREIT figures of merit in a D-Bar setting, as the GREIT figures of merit assume that the reconstruction background is centered around 0, which is the case for the majority of time-differential EIT algorithms, such as GREIT, Gauss–Newton, or NOSER. The D-Bar reconstructions are centered around 1; thus, impedance inclusions with a conductivity greater than the background are greater than one, and conductivities smaller than the background are smaller than one, but still above zero. We simply subtract 1 from the reconstructed impedance values to apply the GREIT figures of merit, as this allows us to use the standard GREIT figures of merit. The figures of merit are as follows:
Firstly, the quarter-amplitude set x q is calculated, which is defined as
x ^ q = { 1 [ x ] i 1 4 m a x ( x ) 0 o t h e r w i s e ,
where [ x ] i is each pixel in the reconstructed image. According to Adler [36], the value of the threshold is somewhat arbitrary and has little effect on the results. Amplitude response (AR) is the summation of all pixels in the quarter-amplitude set divided by the sum of all pixels in the reference image. It is calculated as follows:
AR = x ^ q x 0 ,
where x 0 is the impedance target pixels of the reference image. The AR should be close to 1.
The position error (PE) measures the distance from the center of gravity ( c o g ) of the true conductivity to that of the reconstructed conductivity. Here, we differ slightly from the original definition as our domain is not symmetric; however, the general idea is still the same. We define the PE as
PE = ( c o g x o c o g ^ x ) 2 + ( c o g y o c o g ^ y ) 2 ,
where c o g i o i [ 1 , 2 ] is the i t h component of the center of gravity of the true conductivity distribution and c o g ^ i i [ 1 , 2 ] is the i t h component of the center of gravity of the reconstructed quarter-amplitude-set conductivity. The PE should be as small as possible.
The reconstruction resolution (RES) measures the fraction of the area of the quarter-amplitude set to the area of the entire medium. It is calculated as
RES = A q A 0 ,
where A q is defined as the area of the quarter-amplitude set and A 0 is defined as the area of the entire medium. The RES should be low, as this allows one to distinguish impedance enclosures better.
The shape deformation (SD) measures the fraction of the quarter-amplitude set which does not fit into an equal-sized circle area. It is defined as
SD = j C [ x ^ q ] j j [ x ^ q ] j
The SD should be low, as the high values indicate inhomogeneous reconstructions.
Ringing (RNG) measures overshooting of the reconstruction. It is defined as the area with a negative sign outside the reconstruction. It is the ratio of negative amplitude outside the equivalent circle to the area of amplitude within the equivalent circle. It is defined as
RNG = j C & x ^ < 0 [ x ^ ] j j C [ x ^ ] j
RNG should be low.

4. Results

Firstly, we gave an example of an EIT reconstruction in our squared domain. We used the D-Bar algorithm described in the methods to obtain the reconstructions. The reconstruction can be observed in Figure 8, in the bottom left. The EIT reconstruction of the target is bigger than the ground truth. Throughout the whole paper, the truncation radius was chosen empirically, in order tobetter see the effects. To also show the performance of the reconstruction on two targets, we added, in Figure 9, a plot of two targets with the settings. It showed the same characteristic of better resolution along the x-axis. The values were in the range typically reported by the literature [13,22].
As a second reconstruction example, we tested an EIT reconstruction with R 1 = 3 , but added noise of 0.75 % . The ground truth, the EIT reconstruction, and the combination with radar data are shown in Figure 10. Compared to Figure 8, the EIT reconstruction is more blurred, while the combined reconstruction is thinner along the x-axis.
As a third example, we gave two reconstructions with a priori information and two different estimations of the conductivity in the a priori data. The results are shown in Figure 11. The underestimation (scale factor of 0.001 ) resulted in a reconstruction similar to that of Figure 8, while the overestimation is similar to the combined reconstruction in Figure 10.
Regarding the analysis with the GREIT figures of merit, we tested the EIT and radar reconstruction for an EIT truncation radius of R 1 = 3 without noise in Figure 12. AR became lower with an increasing radar truncation radius R 2 , and reached a plateau for R 2 = 4.75 . It then decreased again for R 2 = 5.5 . The AR with radar information was lower than without. PE followed another pattern: it decreased until R 2 = 4.25 was reached and increased after that. The trend for RES followed the same behavior as AR. It first decreased until R 2 = 4.75 . After that, it reached a plateau, but the value was still lower than without radar information. It increased again at R 2 = 5.5 . SD decreased with increasing truncation radius R 2 . Instead of a plateau after R 2 = 4.7 , the value increased again and surpassed the pure EIT reconstruction. RNG behaved differently; instead of decreasing, it increased steadily until it reached a plateau at R 2 = 4 ; after that, it decreased again.
Similar results occurred when adding the noise to the EIT measurements, which can be seen in Figure 13. Thus, the benefits occur also when noise is present.
The graphs of GREIT figures of merit are qualitatively similar to Figure 11; however, when 0.75 % noise is added, increasing R 2 at first does not give additional spatial information, as all the performance metrics stay in a plateau until around R 2 = 4 . Subsequently, the performance metrics have the same behavior as in the no-noise case.
We next analyzed the behavior of wrong conductivity estimations in the radar a priori information. While the conductivity estimate of the a priori radar data for the previous plots was taken from the literature, we altered the conductivity estimate for the radar. The estimated conductivity from the radar measurements were scaled by γ : = { 0.001 , 0.01 , 0.1 , 1 , 10 } . R 1 was set to 3 and R 2 = 4 , while no noise was used. The results can be seen in Figure 14. AR, PE, RES, and SD have the same qualitative behavior. They are on a plateau from 0.001 to 0.01 and start to decrease for λ from 0.01 to 1. From 1 to 100, the values are building a plateau and increasing only slightly. For practical use, Alsaker et al. developed a nonlinear estimation algorithm for the estimation of an a priori conductivity [37].
When fixing R 1 and increasing R 2 , we could observe that, at first, the results tended to be better, but after R 2 = 4.5 , the figures of merit tended to become worse again in some figures and a little bit better in others. At first, this behavior seems strange; however, the radar information is only available along the x-axis. Thus, the more radar information is added, the more the reconstruction tends to be narrowed down along the x-axis, but, at the same time, is expanded along the y-axis. This results in an increase in the PE as the center of gravity of the quarter-amplitude set is drawn to the center of the y-axis. The same holds true for the SD. The behavior of the AR and RES, however, will become better as more radar information is added. This can still be explained. As the radar information contains solely values from the k-space with R 1 < | k | , only the fine-detail components are enteringinto the reconstructed image. Thus, the area of the reconstruction tends to decrease, which results in an increase in the aforementioned metrics; however, the best value for R 2 in our setup lies around 4.0 .
The conductivity for the calculation of the scattering transform needs to be estimated because our radar setup only provides spatial information. In the previous case, we took the exact value of simulated conductivity. This, however, is not a valid estimation for the real-world case. We simulated different conductivity estimates to check the stability. Since the simulated enclosure has a higher value than the background conductivity, overestimating the conductivity tends to lead to slightly better results, which are barely observable in the reconstruction (see Figure 11). However, the improvements are very minor. This might no longer be the case for more complicated geometries. Underestimation of the conductivity, however, leads to a strong decrease in reconstruction quality, as low-frequency spatial information from EIT and high-frequency spatial information from radar in the k-space compete against one another during the reconstruction.

5. Conclusions and Outlook

We showed that radar information along the x-axis can improve EIT reconstructions up to a certain point. The resolution increases along the x-axis; this can be seen especially in Figure 9, as the two targets are better distinguishable. This is due to the fact that the EIT data need to be truncated in the k a -space, while the radar data are less noisy. Thus, the added information in the k-space introduces more high-frequency components into the reconstruction, which leads to more detail in the resulting images.
Adding too much radar information in this setting can result in a worse performance of the reconstruction algorithm. When trying to reconstruct a single object with conductivity higher than the background, an overestimation of the conductivity tends to be more forgiving than underestimating the conductivity.
The shape of the truncation in the k-space needs to be investigated more in future work. Shapes other than a circle may lead to even better results, and may potentially be able to overcome some disadvantages when only additional information along the x-axis is used. However, it is not clear what this shape might look like, as the k-space representation has no humanly interpretable link to the real domain.
An improvement of the shape correction term used for the DN-map may even further improve the results.
A more sophisticated construction of the prior might also lead to better results; for this paper, however, we followed the steps for the EIT reconstruction as it is easier to follow and does not change the findings.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/a16010043/s1.

Author Contributions

Conceptualization, J.R., J.M., and C.N.; methodology, J.R. and C.N.; software, J.R. and D.H.N.; validation, J.R., J.M., and C.N.; formal analysis, J.R. and D.H.N.; investigation, J.R. and S.L.; resources, S.L. and J.M.; data curation, J.M. and D.H.N.; writing—original draft preparation, J.R.; writing—review and editing, J.R. and C.N.; visualization, J.R.; supervision, C.N. and S.L.; project administration, S.L.; funding acquisition, J.M. and S.L. All authors have read and agreed to the published version of the manuscript.

Funding

The authors gratefully acknowledge financial support provided by the project InDiThera (13GW0361D, 13GW0361E) through the Federal Ministry of Education and Research.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
EITElectrical impedance tomography
DN-mapDirichlet-to-Neumann map
NDNeumman-to-Dirichlet
RMSDRoot-mean-square deviation
ARAmplitude response
PEPosition error
RESResolution
SDShape deformation
RNGRinging
CGOComplex geometric optics

Appendix A

Appendix A.1. D-Bar Algorithm

The D-Bar method uses complex geometrical optics (CGO) solutions to solve the conductivity σ inside Ω , using the DN-map as an input. As a first step, the Laplace Equation (1) is transformed to the time-independent Schrödinger equation, which is achieved through the change of variables. With q ( z ) = σ 1 / 2 Δ σ 1 / 2 and v = σ 1 / 2 u , Equation (1) becomes
( Δ + q ( z ) ) v = 0 .
Using the assumption that σ is constant near the boundary Ω , the now more abstract potential q ( z ) is extended from the domain Ω to the whole complex plane with q ( z ) = 0 outside Ω . This enables the integration along the boundary to be simpler. The CGO solutions can now be used to solve q ( z ) . After introducing a parameter k, we define the solution by Ψ ( k , z ) —both k and Ψ are complex-valued—which satisfies
( Δ + q ( z ) ) Ψ ( k , z ) = 0 .
These CGO solutions were used by L. D. Faddeev in the context of quantum scattering [38]. Subsequently, Astala and Päivärintas used CGO solutions for their solution of the Calderón problem (to determine the conductivity given the DN map) [27]. In order to solve for Ψ of Equation (A2), we define a function μ ( k , z ) by
μ ( k , z ) : = e j k z Ψ ( k , z ) .
It is defined this way because Ψ is asymptotic to e i k z and it also satisfies the ¯ (D-Bar) equation, which gave the algorithm its name. It is given in integral form as [28]
μ ( k , z ) = 1 + 1 ( 2 π ) 2 R 2 t ( k ) ( s k ) k ¯ e z μ ( k , z ) ¯ d k 1 d k 2 .
Here, e z = e j ( k z + k ¯ z ¯ ) and t ( k ) is the so-called scattering transform. We recommend reading [28] for further details on the mathematical background. The scattering transform can be regarded as a nonlinear Fourier transformation, i.e., it transforms inputs to a complex frequency-like domain with nonlinear waves. It is defined as
t ( k ) = Ω e i k ¯ z ¯ ( Λ σ Λ 1 ) Ψ ( k , z ) d s ,
where Λ 1 is the DN map for the case of conductivity 1 inside the domain Ω and d s is the arc length of the domain boundary. When the D-Bar equation is solved for μ ( k , z ) , one can calculate the conductivity inside the domain by
lim k 0 μ ( k , z ) 2 = σ ( z ) , z Ω
Note that (A6) is dependent on the z-position. Therefore, the D-Bar algorithm can solve the conductivity for specific points inside the domain. The density of reconstruction points can be selected to balance the computation speed and the potential spatial resolution. Another option is to just reconstruct points in a certain region of interest. We recommend reading [28] for a thorough explanation of the mathematical background.
The steps described above have a disadvantage in the case of real-world measurements. When the data are noisy, the reconstruction results in distortions. One can imagine the space of the scattering transform as a complex frequency-like domain, where | k | can be interpreted as a kind of frequency. High-frequency (larger | k | ) components provide greater detail in terms of spatial resolution, for example, for sharp edges, but are simultaneously more prone to noise. Thus, the scattering transformation t ( k ) is not evaluated in the whole complex plane, but, rather, on a truncated disk of radius | k | < = R 1 around the origin. The scattering transformation outside R 1 < | k | is set to 0.

Appendix A.2. Sample of K

Figure A1. This figure shows a simplified example of the content of K . In this, M = 4 (which implies m = 2 ) and a truncation radius R = 4 was chosen. This grid is used through steps 2 and 3 of the D-Bar algorithm.
Figure A1. This figure shows a simplified example of the content of K . In this, M = 4 (which implies m = 2 ) and a truncation radius R = 4 was chosen. This grid is used through steps 2 and 3 of the D-Bar algorithm.
Algorithms 16 00043 g0a1

References

  1. Zhu, Q.; Lionheart, W.; Lidgey, F.; McLeod, C.; Paulson, K.; Pidcock, M. An adaptive current tomograph using voltage sources. IEEE Trans. Biomed. Eng. 1993, 40, 163–168. [Google Scholar] [CrossRef] [PubMed]
  2. Hentze, B.; Muders, T.; Luepschen, H.; Maripuu, E.; Hedenstierna, G.; Putensen, C.; Walter, M.; Leonhardt, S. Regional lung ventilation and perfusion by electrical impedance tomography compared to single-photon emission computed tomography. Physiol. Meas. 2018, 39, 065004. [Google Scholar] [CrossRef] [PubMed]
  3. Leonhardt, S.; Cordes, A.; Plewa, H.; Pikkemaat, R.; Soljanik, I.; Moehring, K.; Gerner, H.J.; Rupp, R. Electric impedance tomography for monitoring volume and size of the urinary bladder. Biomed. Eng./Biomed. Tech. 2011, 56, 301–307. [Google Scholar] [CrossRef] [PubMed]
  4. Meier, T.; Luepschen, H.; Karsten, J.; Leibecke, T.; Großherr, M.; Gehring, H.; Leonhardt, S. Assessment of regional lung recruitment and derecruitment during a PEEP trial based on electrical impedance tomography. Intensive Care Med. 2008, 34, 543–550. [Google Scholar] [CrossRef] [PubMed]
  5. Abascal, J.F.P.; Arridge, S.R.; Atkinson, D.; Horesh, R.; Fabrizi, L.; De Lucia, M.; Horesh, L.; Bayford, R.H.; Holder, D.S. Use of anisotropic modelling in electrical impedance tomography; description of method and preliminary assessment of utility in imaging brain function in the adult human head. Neuroimage 2008, 43, 258–268. [Google Scholar] [CrossRef]
  6. Hong, S.; Lee, K.; Ha, U.; Kim, H.; Lee, Y.; Kim, Y.; Yoo, H.J. A 4.9 mΩ-sensitivity mobile electrical impedance tomography IC for early breast-cancer detection system. IEEE J. Solid-State Circuits 2014, 50, 245–257. [Google Scholar] [CrossRef]
  7. Kruger, M.; Poolla, K.; Spanos, C.J. A class of impedance tomography based sensors for semiconductor manufacturing. In Proceedings of the 2004 American Control Conference, Boston, MA, USA, 30 June–2 July 2004; Volume 3, pp. 2178–2183. [Google Scholar]
  8. Hallaji, M.; Seppänen, A.; Pour-Ghaz, M. Electrical impedance tomography-based sensing skin for quantitative imaging of damage in concrete. Smart Mater. Struct. 2014, 23, 085001. [Google Scholar] [CrossRef]
  9. Vauhkonen, M.; Vadasz, D.; Karjalainen, P.A.; Somersalo, E.; Kaipio, J.P. Tikhonov regularization and prior information in electrical impedance tomography. IEEE Trans. Med. Imaging 1998, 17, 285–293. [Google Scholar] [CrossRef]
  10. Borsic, A.; Graham, B.M.; Adler, A.; Lionheart, W.R. Total Variation Regularization in Electrical Impedance Tomography; University of Manchester: Manchester, UK, 2007. [Google Scholar]
  11. Kaipio, J.P.; Kolehmainen, V.; Somersalo, E.; Vauhkonen, M. Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography. Inverse Probl. 2000, 16, 1487. [Google Scholar] [CrossRef] [Green Version]
  12. Ahmad, S.; Strauss, T.; Kupis, S.; Khan, T. Comparison of statistical inversion with iteratively regularized Gauss Newton method for image reconstruction in electrical impedance tomography. Appl. Math. Comput. 2019, 358, 436–448. [Google Scholar] [CrossRef]
  13. Isaacson, D.; Mueller, J.L.; Newell, J.C.; Siltanen, S. Reconstructions of chest phantoms by the D-bar method for electrical impedance tomography. IEEE Trans. Med. Imaging 2004, 23, 821–828. [Google Scholar] [CrossRef]
  14. Calderón, A.P. On an inverse boundary value problem. Comput. Appl. Math 2006, 25, 133–138. [Google Scholar] [CrossRef] [Green Version]
  15. Nachman, A.I. Global uniqueness for a two-dimensional inverse boundary value problem. Ann. Math. 1996, 143, 71–96. [Google Scholar] [CrossRef]
  16. Siltanen, S.; Mueller, J.; Isaacson, D. An implementation of the reconstruction algorithm of A Nachman for the 2D inverse conductivity problem. Inverse Probl. 2000, 16, 681. [Google Scholar] [CrossRef]
  17. Murphy, E.K.; Mueller, J.L. Effect of domain shape modeling and measurement errors on the 2-D D-bar method for EIT. IEEE Trans. Med. Imaging 2009, 28, 1576–1584. [Google Scholar] [CrossRef]
  18. Sun, B.; Yue, S.; Hao, Z.; Cui, Z.; Wang, H. An improved Tikhonov regularization method for lung cancer monitoring using electrical impedance tomography. IEEE Sens. J. 2019, 19, 3049–3057. [Google Scholar] [CrossRef]
  19. Adler, A.; Lionheart, W.R. Uses and abuses of EIDORS: An extensible software base for EIT. Physiol. Meas. 2006, 27, S25. [Google Scholar] [CrossRef] [Green Version]
  20. Soleimani, M. Electrical impedance tomography imaging using a priori ultrasound data. Biomed. Eng. Online 2006, 5, 8. [Google Scholar] [CrossRef]
  21. Alsaker, M.; Mueller, J.L. A D-bar algorithm with a priori information for 2-dimensional electrical impedance tomography. SIAM J. Imaging Sci. 2016, 9, 1619–1654. [Google Scholar] [CrossRef]
  22. Hamilton, S.J.; Mueller, J.L.; Alsaker, M. Incorporating a Spatial Prior into Nonlinear D-Bar EIT Imaging for Complex Admittivities. IEEE Trans. Med. Imaging 2017, 36, 457–466. [Google Scholar] [CrossRef]
  23. Adler, A.; Gaggero, P.; Maimaitijiang, Y. Distinguishability in EIT using a hypothesis-testing model. J. Phys. Conf. Ser. 2010, 224, 012056. [Google Scholar] [CrossRef] [Green Version]
  24. Prati, M.V.; Moll, J.; Kexel, C.; Nguyen, D.H.; Santra, A.; Aliverti, A.; Krozer, V.; Issakov, V. Breast cancer imaging using a 24 GHz ultra-wideband MIMO FMCW radar: System considerations and first imaging results. In Proceedings of the 2020 14th European Conference on Antennas and Propagation (EuCAP), Copenhagen, Denmark, 15–20 March 2020; pp. 1–5. [Google Scholar]
  25. Flores-Tapia, D.; Pistorius, S. Electrical impedance tomography reconstruction using a monotonicity approach based on a priori knowledge. In Proceedings of the 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology, Buenos Aires, Argentina, 31 August–4 September 2010; pp. 4996–4999. [Google Scholar]
  26. Cheney, M.; Isaacson, D. Issues in electrical impedance imaging. IEEE Comput. Sci. Eng. 1995, 2, 53–62. [Google Scholar] [CrossRef]
  27. Astala, K.; Päivärinta, L. Calderón’s inverse conductivity problem in the plane. Ann. Math. 2006, 163, 265–299. [Google Scholar] [CrossRef] [Green Version]
  28. Mueller, J.L.; Siltanen, S. Linear and Nonlinear Inverse Problems with Practical Applications; SIAM: Philadelphia, PA, USA, 2012. [Google Scholar]
  29. Herrera, C.N.; Vallejo, M.F.; Mueller, J.L.; Lima, R.G. Direct 2-D reconstructions of conductivity and permittivity from EIT data on a human chest. IEEE Trans. Med. Imaging 2014, 34, 267–274. [Google Scholar] [CrossRef] [PubMed]
  30. Dodd, M.; Mueller, J.L. A real-time D-bar algorithm for 2-D electrical impedance tomography data. Inverse Probl. Imaging 2014, 8, 1013. [Google Scholar] [CrossRef]
  31. Saad, Y.; Schultz, M.H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 1986, 7, 856–869. [Google Scholar] [CrossRef] [Green Version]
  32. Gabriel, S.; Lau, R.; Gabriel, C. The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues. Phys. Med. Biol. 1996, 41, 2271. [Google Scholar] [CrossRef] [Green Version]
  33. Surowiec, A.J.; Stuchly, S.S.; Barr, J.R.; Swarup, A. Dielectric properties of breast carcinoma and the surrounding tissues. IEEE Trans. Biomed. Eng. 1988, 35, 257–263. [Google Scholar] [CrossRef]
  34. Martellosio, A.; Pasian, M.; Bozzi, M.; Perregrini, L.; Mazzanti, A.; Svelto, F.; Summers, P.E.; Renne, G.; Preda, L.; Bellomi, M. Dielectric properties characterization from 0.5 to 50 GHz of breast cancer tissues. IEEE Trans. Microw. Theory Tech. 2016, 65, 998–1011. [Google Scholar] [CrossRef]
  35. Mueller, J.L.; Siltanen, S. The d-bar method for electrical impedance tomography—Demystified. Inverse Probl. 2020, 36, 093001. [Google Scholar] [CrossRef]
  36. Adler, A.; Arnold, J.H.; Bayford, R.; Borsic, A.; Brown, B.; Dixon, P.; Faes, T.J.; Frerichs, I.; Gagnon, H.; Gärber, Y.; et al. GREIT: A unified approach to 2D linear EIT reconstruction of lung images. Physiol. Meas. 2009, 30, S35. [Google Scholar] [CrossRef]
  37. Alsaker, M.; Mueller, J.L. Use of an optimized spatial prior in D-bar reconstructions of EIT tank data. Inverse Probl. Imaging 2018, 12, 883. [Google Scholar] [CrossRef]
  38. Faddeev, L.D. Increasing solutions of the Schrödinger equation. Sov. Phys. Dokl. 1966, 10, 1033. [Google Scholar]
Figure 1. Visualization of one stimulation pattern. Ω surrounded by the thick black circle denotes the domain of interest. The red line indicates the current depending on the angle Θ . The dashed gray line denotes a current amplitude of zero. Note that this is not the domain used in this paper, but a standard circular domain for better visualization of the concept.
Figure 1. Visualization of one stimulation pattern. Ω surrounded by the thick black circle denotes the domain of interest. The red line indicates the current depending on the angle Θ . The dashed gray line denotes a current amplitude of zero. Note that this is not the domain used in this paper, but a standard circular domain for better visualization of the concept.
Algorithms 16 00043 g001
Figure 2. Visualization of the finite element method model with the electrode position in cyan and domain boundary in blue. r ( Θ 1 ) and r ( Θ 2 ) are the positional angles of the first and second electrodes. The angle difference Δ Θ 2 , 1 is needed for the calculation of the derivative in Equation (6).
Figure 2. Visualization of the finite element method model with the electrode position in cyan and domain boundary in blue. r ( Θ 1 ) and r ( Θ 2 ) are the positional angles of the first and second electrodes. The angle difference Δ Θ 2 , 1 is needed for the calculation of the derivative in Equation (6).
Algorithms 16 00043 g002
Figure 3. Example of the k-space. At the top, the conductivity enclosure in our rectangular domain is given. In the bottom left part, the real part of the k-space is given, while the imaginary part is displayed on the right. The truncation radius is R 1 = 3.5 . It is important to note that the size of the actual domain has no direct relation to the more abstract truncation radius, as the truncation radius is placed in the k-space.
Figure 3. Example of the k-space. At the top, the conductivity enclosure in our rectangular domain is given. In the bottom left part, the real part of the k-space is given, while the imaginary part is displayed on the right. The truncation radius is R 1 = 3.5 . It is important to note that the size of the actual domain has no direct relation to the more abstract truncation radius, as the truncation radius is placed in the k-space.
Algorithms 16 00043 g003
Figure 4. Example of the sample strategy of K for p = 2 which results in a 4 × 4 grid (note that this is just for illustration purposes and would be too low for computations). Every quadrant of the k-space has the same number of points. The numbers above the points are the indices of the matrix K as implemented in Matlab®.
Figure 4. Example of the sample strategy of K for p = 2 which results in a 4 × 4 grid (note that this is just for illustration purposes and would be too low for computations). Every quadrant of the k-space has the same number of points. The numbers above the points are the indices of the matrix K as implemented in Matlab®.
Algorithms 16 00043 g004
Figure 5. Structure of the used FEM model. On the bottom and top are the dual-use EIT electrodes/radar antennas. On the right and left are pure EIT electrodes. The skin is colored in purple. The rest of the breast is modeled with fat, which is represented by the pink-colored areas. The size of the model is 100 mm × 100 mm. This is the same as in Prati et al. [24].
Figure 5. Structure of the used FEM model. On the bottom and top are the dual-use EIT electrodes/radar antennas. On the right and left are pure EIT electrodes. The skin is colored in purple. The rest of the breast is modeled with fat, which is represented by the pink-colored areas. The size of the model is 100 mm × 100 mm. This is the same as in Prati et al. [24].
Algorithms 16 00043 g005
Figure 6. Reconstruction of two tumors using the reconstruction algorithm described. The RMSD curve can be seen. The gray dashed line is the envelope that aids in the calculation of the estimated tumor position, which is the red line—adapted from [24]. The position of the tumor along the y-axis is 27.5 mm and it has a diameter of 8 mm.
Figure 6. Reconstruction of two tumors using the reconstruction algorithm described. The RMSD curve can be seen. The gray dashed line is the envelope that aids in the calculation of the estimated tumor position, which is the red line—adapted from [24]. The position of the tumor along the y-axis is 27.5 mm and it has a diameter of 8 mm.
Algorithms 16 00043 g006
Figure 7. Visualization of the merging in the k-space. The purple area represents the k-space that is filled with t radar ( k ) , while the red area is filled with t E I T ( k ) . Values outside | k | < = R 2 are zero.
Figure 7. Visualization of the merging in the k-space. The purple area represents the k-space that is filled with t radar ( k ) , while the red area is filled with t E I T ( k ) . Values outside | k | < = R 2 are zero.
Algorithms 16 00043 g007
Figure 8. Top: Ground truth of the conductivity enclosure. Bottom left: EIT reconstruction with no noise added and R 1 = 3 . Note: There is no a priori information in the image. Bottom right: Combined reconstruction with a priori information with a truncation radius R 2 = 4 .
Figure 8. Top: Ground truth of the conductivity enclosure. Bottom left: EIT reconstruction with no noise added and R 1 = 3 . Note: There is no a priori information in the image. Bottom right: Combined reconstruction with a priori information with a truncation radius R 2 = 4 .
Algorithms 16 00043 g008
Figure 9. Top: Ground truth of the two conductivity enclosures. Bottom left: EIT reconstruction with no noise added and R 1 = 3 . Note: There is no a priori information in the image. Bottom right: Combined reconstruction with a priori information with a truncation radius R 2 = 4 .
Figure 9. Top: Ground truth of the two conductivity enclosures. Bottom left: EIT reconstruction with no noise added and R 1 = 3 . Note: There is no a priori information in the image. Bottom right: Combined reconstruction with a priori information with a truncation radius R 2 = 4 .
Algorithms 16 00043 g009
Figure 10. Top: Ground truth of the conductivity enclosure. Bottom left: EIT reconstruction with 0.75 % noise added and R 1 = 3 . Bottom right: Combination of noisy EIT data ( 0.75 % noise) and radar data. R 1 = 3 and R 2 = 4 . The noise was added to the voltages of the electrodes in percentage of the peak-to-peak value.
Figure 10. Top: Ground truth of the conductivity enclosure. Bottom left: EIT reconstruction with 0.75 % noise added and R 1 = 3 . Bottom right: Combination of noisy EIT data ( 0.75 % noise) and radar data. R 1 = 3 and R 2 = 4 . The noise was added to the voltages of the electrodes in percentage of the peak-to-peak value.
Algorithms 16 00043 g010
Figure 11. Visualization of under- and overestimation of the gourd truth conductivity. Top: Ground truth of the conductivity enclosure. Bottom left: Combination reconstruction with R 1 = 3 and R 2 = 4 . The original conductivity of the enclosure is underestimated by a factor of 0.001 . Bottom right: Combination reconstruction with R 1 = 3 and R 2 = 4 . The original conductivity of the enclosure is overestimated by a factor of 100.
Figure 11. Visualization of under- and overestimation of the gourd truth conductivity. Top: Ground truth of the conductivity enclosure. Bottom left: Combination reconstruction with R 1 = 3 and R 2 = 4 . The original conductivity of the enclosure is underestimated by a factor of 0.001 . Bottom right: Combination reconstruction with R 1 = 3 and R 2 = 4 . The original conductivity of the enclosure is overestimated by a factor of 100.
Algorithms 16 00043 g011
Figure 12. GREIT figures of merit for the combination of EIT and radar data, for R 1 = 3 without noise. The orange line represents the scores of the EIT-only data. R 2 was increased from 3.25 to 5.5 in steps of 0.25 .
Figure 12. GREIT figures of merit for the combination of EIT and radar data, for R 1 = 3 without noise. The orange line represents the scores of the EIT-only data. R 2 was increased from 3.25 to 5.5 in steps of 0.25 .
Algorithms 16 00043 g012
Figure 13. GREIT figures of merit for the combination of EIT and radar data, for R 1 = 3 with 0.75 % noise. The orange line represents the scores of the EIT-only data. R 2 was increased from 3.25 to 5.5 in steps of 0.25 .
Figure 13. GREIT figures of merit for the combination of EIT and radar data, for R 1 = 3 with 0.75 % noise. The orange line represents the scores of the EIT-only data. R 2 was increased from 3.25 to 5.5 in steps of 0.25 .
Algorithms 16 00043 g013
Figure 14. GREIT figures of merit for the combination of EIT and radar data, for R 1 = 3 without noise and R 2 = 4 . The radar estimation conductivity was changed and scaled by γ : = { 0.001 , 0.01 , 0.1 , 1 , 10 } . The x-axis is scaled logarithmically.
Figure 14. GREIT figures of merit for the combination of EIT and radar data, for R 1 = 3 without noise and R 2 = 4 . The radar estimation conductivity was changed and scaled by γ : = { 0.001 , 0.01 , 0.1 , 1 , 10 } . The x-axis is scaled logarithmically.
Algorithms 16 00043 g014
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Rixen, J.; Leonhardt, S.; Moll, J.; Nguyen, D.H.; Ngo, C. The D-Bar Algorithm Fusing Electrical Impedance Tomography with A Priori Radar Data: A Hands-On Analysis. Algorithms 2023, 16, 43. https://doi.org/10.3390/a16010043

AMA Style

Rixen J, Leonhardt S, Moll J, Nguyen DH, Ngo C. The D-Bar Algorithm Fusing Electrical Impedance Tomography with A Priori Radar Data: A Hands-On Analysis. Algorithms. 2023; 16(1):43. https://doi.org/10.3390/a16010043

Chicago/Turabian Style

Rixen, Jöran, Steffen Leonhardt, Jochen Moll, Duy Hai Nguyen, and Chuong Ngo. 2023. "The D-Bar Algorithm Fusing Electrical Impedance Tomography with A Priori Radar Data: A Hands-On Analysis" Algorithms 16, no. 1: 43. https://doi.org/10.3390/a16010043

APA Style

Rixen, J., Leonhardt, S., Moll, J., Nguyen, D. H., & Ngo, C. (2023). The D-Bar Algorithm Fusing Electrical Impedance Tomography with A Priori Radar Data: A Hands-On Analysis. Algorithms, 16(1), 43. https://doi.org/10.3390/a16010043

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop