Analytical and numerical models of the magnetoacoustic tomography with magnetic induction

Marcin Ziolkowski (Faculty of Electrical Engineering, West Pomeranian University of Technology, Szczecin, Poland)
Stanislaw Gratkowski (Faculty of Electrical Engineering, West Pomeranian University of Technology, Szczecin, Poland)
Adam Ryszard Zywica (Faculty of Electrical Engineering, West Pomeranian University of Technology, Szczecin, Poland)

Abstract

Purpose

Electrical properties of biological tissues are known to be sensitive to physiological and pathological conditions of living organisms. For instance, human breast cancer or liver tumor cells have a significantly higher electrical conductivity than a healthy tissue. The paper aims to the new recently developed magnetoacoustic tomography with magnetic induction (MAT-MI) which can be deployed for electrical conductivity imaging of low-conductivity objects. Solving a test problem by using an analytical method is a useful exercise to check the validity of the more complex numerical finite element models. Such test problems are discussed in Chapter 3. The detailed analysis of an electromagnetic induction in low-conductivity objects is very important for the next steps in the tomographic process of image reconstruction. Finally, the image reconstruction examples for object’s complex shapes’ have been analyzed. The Lorentz force divergence reconstruction has been achieved with the help of time reversal algorithm.

Design/methodology/approach

In given arrangements the magnetic field and eddy current vectors satisfy the Maxwell partial differential equations. Applying the separation of variables method analytical solutions are obtained for an infinitely long conducting cylindrical segment in transient magnetic field. A special case for such a configuration is an infinitely long cylinder with longitudinal crack. The analytical solutions are compared with those obtained by using numerical procedures. For complex shapes of the object, the MAT-MI images have been calculated with the help of the finite element method and time reversal algorithm.

Findings

The finite element model developed for a MAT-MI forward problem has been validated by analytical formulas. Based on such a confirmation, the MAT-MI complex model has been defined and solved. The conditions allowing successful MAT-MI image reconstruction have been provided taking into account different conductivity distribution. For given object’s parameters, the minimum number of measuring points allowing successful reconstruction has been determined.

Originality/value

A simple test example has been proposed for MAT-MI forward problem. Analytical closed-form solutions have been used to check the validity of the made in-house finite element software. More complex forward and inverse problems have been solved using the software.

Keywords

Citation

Ziolkowski, M., Gratkowski, S. and Zywica, A. (2018), "Analytical and numerical models of the magnetoacoustic tomography with magnetic induction", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. 37 No. 2, pp. 538-548. https://doi.org/10.1108/COMPEL-12-2016-0530

Download as .RIS

Publisher

:

Emerald Publishing Limited

Copyright © 2018, Emerald Publishing Limited


1. Introduction

Many various imaging techniques are used in contemporary biomedical research and clinical diagnosis. A new recently developed method is called magnetoacoustic tomography with magnetic induction (MAT-MI). Historically, MAT-MI was first proposed by Xu and He (2005). This tomographic technique is especially dedicated for a non-invasive imaging of low-conductivity objects (e.g. biological tissues) and bases on coupling ultrasound and magnetism into one unified hybrid imaging system. In such a method, an object is placed in static and time-varying (pulsed) magnetic fields. The interaction between induced eddy currents and static magnetic field leads to generate the Lorentz force. The divergence of the Lorentz force acts as an ultrasonic waves’ source, and the transmitted acoustic signals are acquired by set of piezoelectric transducers placed around the object. Finally, the conductivity map can be reconstructed with the help of various back projection methods. Combining the good contrast of electrical impedance tomography and the high spatial resolution of ultrasonography, MAT-MI seems to be promising technique in medical imaging diagnosis (Li et al., 2016; Ma and He, 2008; Sun et al., 2015).

MAT-MI method involves two main parts. The first, so-called forward problem comes down to determine the ultrasonic wave as a change of acoustic pressure originated in the object of low-conductivity. The latter, called inverse problem, consists of two steps, i.e. reconstructing the distribution of the Lorentz force divergence according to the ultrasonic signals collected by piezoelectric transducers, and finally the imaging of the electrical conductivity distribution (Ammari et al., 2015; Żywica, 2016).

The detailed analysis of an electromagnetic induction in low-conductivity objects is very important for the next steps in the tomographic process of image reconstruction. In this paper, we discuss, very shortly, the MAT-MI principle. Next, the first step of the MAT-MI forward problem, i.e. analysis of the transient eddy currents is performed. Numerical simulations of the transient eddy current phenomenon are compared with the analytical solutions. Solving a test problem by using an analytical method is a useful exercise to check the validity of the more complex numerical models. Finally, the image reconstruction examples for object’s complex shape are analyzed. The Lorentz force divergence reconstruction is achieved with the help of time reversal algorithm. For a given object’s parameters, the minimum number of measuring points allowing successful reconstruction is determined.

2. Magnetoacoustic tomography with magnetic induction governing theory

The illustration of MAT-MI concept has been shown in Figure 1. Let us consider an object with two layers of different electrical conductivities σ (inner layer, which can be a muscle or brain tissue equivalent) and σ 0 (outer layer, skull or fat equivalent) placed in a static B 0 and time-varying B 1(t) (pulsed) external magnetic fields. The time-varying magnetic field induces eddy currents in the object’s volume. Accordingly, the Lorentz force (result of interaction between a static magnetic field and eddy currents) is generated in the object and finally acoustic vibrations are emitted. Ultrasonic signals acquired by a set of piezoelectric transducers are then used for electrical conductivity image reconstruction process.

The complete MAT-MI forward problem consists of two physical phenomena, i.e. electromagnetic induction in the conductive object and the acoustic wave propagation (result of the interaction between the static magnetic field and the eddy currents). The magnetic stimulation B 1(r , t) (r is the position vector) can be represented by the curl of its magnetic vector potential A (r , t): B (r , t) = ∇ × A (r , t). The conductivity values are small enough that we can ignore the effect of the secondary magnetic field. Additionally, in MAT-MI for driving the stimulating coil ms or µs level current pulses are used and hence the magnetic induction problem can be considered as quasi-static (the displacement currents are neglected). The electric field inside the object and the eddy currents density can be expressed as:

(1)E(r,t)=A(r,t)tΦ(r,t),J(r,t)=σ(r)E(r,t)
where Φ(r , t) is the electric scalar potential.

Taking into account the fact that the magnetic stimulation pulse is short enough, after separating the time-varying magnetic fields for the functions of space and time, the Lorentz force density acting on eddy currents over unit volume can be written as:

(2)F(r,t)=[J(r)×B0(r)]δ(t)
where δ(t) is the Dirac delta function in the time domain.

In the MAT-MI, the divergence of the force resulting from equation (2) leads to mechanical vibrations and generates an acoustic wave governed by the following wave equation:

(3)2p(r,t)1cs22p(r,t)t2=[J(r)×B0(r)]δ(t)
where p(r , t) is the acoustic pressure at spatial point r , cs = (ρ 0·βs)−1/2 is the acoustic speed, ρ 0 is the density of the object at rest and βs is the adiabatic compressibility of the medium.

By using Green’s function, equation (3) can be solved as (Li et al., 2016; Ma and He, 2008):

(4)p(r,t)=14πVr[J(r)×B0(r)]δ(tRcs)R dr
where R = |r r ′|, and the integration is over the sample volume containing the acoustic source.

3. Transient eddy currents in conductive cylindrical segment

Analytic treatment of the eddy-current problems in the cylindrical conductors has been already elaborated quite extensively in the literature (Wwedensky, 1921; Kahn and Spal, 1979). The cited references should only be taken as illustrative examples. Wwedensky developed the solution for the transient magnetic field inside an infinitely long conducting rod surrounded by an infinitely long solenoid, and Khan considered the sinusoidal eddy-current problems in the cylindrical structure possessing the crack.

Let us consider an infinitely long cylinder with the non-regular shape of the cross section. The electrical conductivity and magnetic permeability of the material are σ and µ, respectively (Figure 2). The question is: what is the magnetic field H in the cylindrical segment, if the current i(t) flowing in the winding (which tightly surrounds the segment) is described by the step function, i.e. i(t) = I 11(t)? In consequence, the magnetic field suddenly appears around the cylindrical segment. Mathematically, the phenomenon can be described as a product of the step function and the z–unit vector: H (t) = (nI 1/l)1(t)1 z = H 11(t)1 z, n – number of turns, l – length.

Taking into account the cylindrical coordinate system (r, φ, z), the magnetic field inside the cylinder has only an axial component, H = H(r, φ, t)1 z and satisfies the following partial differential equation:

(5)2Hr2+1rHr+1r22Hφ2=σμHt
with initial and boundary conditions:

H(r, φ, 0) = 0 for r ∈ ⟨(0, a) and φ ∈ ⟨–α, α⟩; H(a, φ, t) = H 1, H(r, ±α, t) = H 1, for t > 0; limtH(r,φ,t)=H1 for t>0 .

The problem can be solved by using the separation of variables method:

(6)H(r,φ,t)=H1+n=1m=0Gm,nJνm(κm,nar)cos(νmφ)exp[κm,n2a2σμt]
where νm=(m  +  0.5)πα , κm , n – successive positive roots of the equation Jνm(κ)=0, m = 0, 1, 2, …, and Gm , n are unknown coefficients which can be found from the initial condition (later, the expression for Gm , n will be given for a special case α = π).

Now, let us consider the situation when α value is equal to π. It means that in the conductive cylinder a longitudinal infinitely thin crack is located (Figure 2, on the right), and the induced currents must bypass such a discontinuity. In this case, it is possible to obtain exact, closed form analytical formula for the induced eddy current J components. The magnetic field inside the cylinder can be determined by following formula:

(7)H(r,φ,t)=H1+n=1m=0Gm,nJm+0.5(κm,nar)cos[(m+0.5)φ]exp(κm,n2a2μσt)
where Jm +0.5(·) denotes Bessel function of the first kind and order m + 0.5 (m = 1, 2, …), respectively; κn (n = 1, 2, …) are the positive roots of the equation Jm +0.5(κ) = 0.

The values of the unknown coefficients Gm , n can be determined by using the expression:

(8)Gm,n=4πH1(1)m2m+10arJm+0.5(κm,nar)dra22Jm0.5(κm,n)Jm+1.5(κm,n)

Taking into account that J   = ∇ × H , the eddy current vector components are:

(9)Jr(r,φ,t)=1rn=1m=0(m+0.5)Gm,nJm+0.5(κm,nar)sin[(m+0.5)φ]exp(κm,n2a2μσt)Jφ(r,φ,t)=1rn=1m=0Gm,n[(m+0.5)Jm+0.5(κm,nar)κm,narJm+1.5(κm,nar)]×  cos[(m+0.5)φ]exp(κm,n2a2μσt)

Figure 3 shows the exemplary eddy current lines for α = π/4 and π, on the left and right, respectively. In Figure 3 on the right, it can be clearly seen that the current strongly turns back in the vicinity of the three characteristic points, i.e. for α ≈ ± π and r ≈ a, and for r ≈ 0.

From the practical point of view, it is essential to compare in quantitative way the results obtained by the formulas (9) with numerical procedures which will be used later for the analysis of complex MAT-MI problem. The calculations have been performed for the electrical conductivity σ = 1 S/m, relative magnetic permeability μr = 1 and cylinder diameter a = 45 mm. In the unit step function, the magnetic field H 1 value has been set to 7.95·105 A/m. The absolute values of the magnetic field H in the cylinder as a function of time t for different r coordinate values have been shown in Figure 4. As expected, in steady state, the magnetic field completely fills the cylinder’s volume. It can be also seen that the analytical and numerical results are in good agreement.

Figure 5 presents the values of the eddy current J as a function of time t for different values of r coordinate. The highest values of the eddy currents can be detected in the vicinity of the cylinder’s circumference. The closer to the middle of the cylinder, the lower values of the current can be observed. It also should be noted that in the steady state eddy currents disappear completely.

Finally, it is convenient to observe the behavior of the eddy current Jφ and Jr components as a function of r coordinate in a given moment of time t for two characteristic values of angle φ. Figure 6 shows the results for the angle φ = 90°. The closer to the circumference of the cylinder the much lower Jr current values can be seen. Approximately for r = 20 mm, the Jφ component changes the sign. It corresponds well to the situation from Figure 6, where for φ = 90° and r ≈ 20 mm, the current strongly turns back.

Another interesting case is for φ very close to the longitudinal crack. Figure 7 presents the eddy current components’ values for φ = 178°.

As expected, the current flows along the crack toward the middle of the conductive cylinder (the lines can be observed in Figure 3), and in consequence the values of the eddy current Jφ are close to zero for all values of r coordinate.

4. Reconstruction of Lorentz force distribution

In this section, we provide examples of acoustic source reconstruction as a Lorentz force distribution for the object of complex shape and complex conductivity distribution. In general, the source term ∇ · (J × B 0) of the wave equation (3) can be solved with the help of well-known time reversal technique. In this case, the main idea of the time reversal method is to reverse in time the recorded acoustic pressure p(r , t) and reemit from the all “measuring points” which are the equivalent of piezoelectric transducers. Next, the resulting vibration will converge back to the point where it was originally emitted. The time reversed signals are described with the following expression (Fink, 1996):

(10)pTR(r,t)=pTR(r,Tt)
where T denotes the duration of this recording.

Let us consider the 2D planar geometry which has been shown in Figure 8. The object to be imaged consists of outer and inner layer placed non-concentrically with low conductivities σ 1 and σ 2, respectively (Beneke et al., 1991). It is assumed that the layers are acoustically homogeneous without any reflections, dispersion and attenuation. The acoustic velocity cs was set to 1490 m/s. The uniform magnetic flux density B 0 value has been set to 1 T. In calculations, it is undertaken that the studies are ideal, and therefore, the recorded signals are considered as a noise-free. The measuring points have been placed on the circumference of the virtual circle with diameter equal to 26 mm.

As there is no analytical solution for the eddy current distribution flowing in the objects of complex shapes, a numerical procedure must be applied. In this case, the MAT-MI forward problem has been performed with the help of finite element method (FEM). Two cases are considered in calculations, i.e. the case “a” (the conductivity σ 2 = 0 S/m, σ 1 > σ 2) and the case “b” (the conductivity σ 2 = 30 S/m, σ 1 < σ 2). Results of the reconstruction for different number of measuring points have been presented in Figures 9 and 10.

It can be observed that in all cases for six measuring points (60° step), some blurring artefacts around and between the object’s region appear (Figures 9 and 10, on the left). However, more measuring points leads to better quality images with less time reversal noise. Significant image reconstruction improvement takes place when the number of measurement points begins by 36 (10° step) and 90 (4° step) for case “a” and case “b”, respectively. Taking into account such requirements, one can obtain good overall contrast pattern with relatively small disturbances located at the objects’ boundaries. It can be also noted that further increment of measuring points’ number does not change the image crucially.

5. Conclusions

In this paper, analytical and numerical models of the magnetoacoustic tomography with magnetic induction have been analyzed and compared. According to the first step of MAT-MI, forward problem analytical formulas for the transient magnetic induction phenomenon have been provided and compared with the finite element method. The obtained results validated successfully numerical models developed with the help of FEM. In case of MAT-MI inverse problem, the object of complex shape and conductivity distribution has been considered. The conditions allowing successful image reconstruction have been provided taking into account different conductivity distribution.

Figures

Illustration of MAT-MI concept

Figure 1.

Illustration of MAT-MI concept

On the left – cross section of the long cylindrical segment with uniform conductivity σ and magnetic permeability µ, on the right – cross section of the long cylinder with longitudinal crack (α = π)

Figure 2.

On the left – cross section of the long cylindrical segment with uniform conductivity σ and magnetic permeability µ, on the right – cross section of the long cylinder with longitudinal crack (α = π)

On the left – exemplary current lines in cylindrical segment for α = π/4; on the right – exemplary current lines for α = π, σ = 1 S/m, t = 10−11 s

Figure 3.

On the left – exemplary current lines in cylindrical segment for α = π/4; on the right – exemplary current lines for α = π, σ = 1 S/m, t = 10−11 s

The absolute values of the magnetic field as a function of time t for φ = 0°, r = 10 and 30 mm

Figure 4.

The absolute values of the magnetic field as a function of time t for φ = 0°, r = 10 and 30 mm

The absolute values of the eddy current J vector as a function of time t for φ = 0°, for different values r

Figure 5.

The absolute values of the eddy current J vector as a function of time t for φ = 0°, for different values r

The values of the eddy current J components for φ = 90° and time t = 10−10 s

Figure 6.

The values of the eddy current J components for φ = 90° and time t = 10−10 s

The values of the eddy current J components for φ = 178° and time t = 10−10 s

Figure 7.

The values of the eddy current J components for φ = 178° and time t = 10−10 s

The object to be imaged with layers of σ
1 and σ
2 and surrounding region σ
0 (distilled water equivalent)

Figure 8.

The object to be imaged with layers of σ 1 and σ 2 and surrounding region σ 0 (distilled water equivalent)

Case “a”: reconstruction of the acoustic source position for six measuring points (on the left) and for 36 measuring points (on the right); σ
0 = 0 S/m, σ
1 = 10 S/m, σ
2 = 0 S/m

Figure 9.

Case “a”: reconstruction of the acoustic source position for six measuring points (on the left) and for 36 measuring points (on the right); σ 0 = 0 S/m, σ 1 = 10 S/m, σ 2 = 0 S/m

Case “b”: reconstruction of the acoustic source position for six measuring points (on the left) and for 90 measuring points (on the right); σ
0 = 0 S/m, σ
1 = 10 S/m, σ
2 = 30 S/m

Figure 10.

Case “b”: reconstruction of the acoustic source position for six measuring points (on the left) and for 90 measuring points (on the right); σ 0 = 0 S/m, σ 1 = 10 S/m, σ 2 = 30 S/m

References

Ammari, H., Boulmier, S. and Millien, P. (2015), “A mathematical and numerical framework for magnetoacoustic tomography with magnetic induction”, Journal of Differential Equations, Vol. 259 No. 10, pp. 5379-5405.

Beneke, R., Neuerburg, J. and Bohndorf, K. (1991), “Muscle cross-section measurement by magnetic resonance imaging”, European Journal of Applied Physiology, Vol. 63 No. 6, pp. 424-429.

Fink, M. (1996), “Time reversal in acoustics”, Contemporary Physics, Vol. 37 No. 2, pp. 95-109.

Kahn, A.H. and Spal, R. (1979), “AC magnetic fields in the vicinity of a crack calculated by analytic and numerical methods”, Proceedings of the DARPA/AFML Review of Progress in Quantitative NDE, July 1978 - September 1979, pp. 65-68.

Li, X., Yu, K. and He, B. (2016), “Magnetoacoustic tomography with magnetic induction (MAT-MI) for imaging electrical conductivity of biological tissue: a tutorial review”, Physics in Medicine & Biology, Vol. 61, pp. 249-270.

Ma, Q. and He, B. (2008), “Magnetoacoustic tomography with magnetic induction: a rigorous theory”, IEEE Transactions on Biomedical Engieneering, Vol. 55 No. 2, pp. 467-479.

Sun, X.D., Wang, X., Zhou, Y.Q., Ma, Q.Y. and Zhang, D. (2015), “Reception pattern influence on magnetoacoustic tomography with magnetic induction”, Chinese Physics B, Vol. 24 No. 1, pp. 014302-014309.

Wwedensky, V.B. (1921), “Concerning the eddy currents generated by a spontaneous change of magnetization”, Annalen Der Physik, Vol. 64, pp. 609-620.

Xu, Y. and He, B. (2005), “Magnetoacoustic tomography with magnetic induction”, Institute of Physics in Medicine and Biology, Vol. 50 No. 21, pp. 5175-5187.

Żywica, A.R. (2016), “Magnetoacoustic tomography with magnetic induction for biological tissue imaging: numerical modelling and simulations”, Archives of Electrical Engineering, Vol. 65 No. 1, pp. 141-150.

Supplementary materials

9781787561342.pdf (42.7 MB)
COMPEL_37_2.pdf (42.8 MB)

Corresponding author

Marcin Ziolkowski is the corresponding author and can be contacted at: marcin.ziolkowski@zut.edu.pl

About the authors

Marcin Ziolkowski received MSc, PhD and DSc (Habilitation) degrees in Electrical Engineering from the Szczecin University of Technology, Szczecin, Poland and West Pomeranian University of Technology, Szczecin, Poland in 2001, 2006 and 2012, respectively. Since 2008, he has been working with the Department of Electrical and Computer Engineering, Electrical Engineering Faculty, West Pomeranian University of Technology, Szczecin, Poland. His main research interests include numerical calculations and visualization of EM fields, inverse problems, electromagnetic field shielding and non-destructive testing of materials.

Stanislaw Gratkowski received MSc, PhD and DSc (Habilitation) degrees in Electrical Engineering from the Szczecin University of Technology, Szczecin, Poland, Gdansk University of Technology, Gdansk, Poland and Electrotechnical Institute, Warsaw, Poland in 1976, 1980 and 1997, respectively. Currently, he is a Professor with the Electrical Engineering Faculty, West Pomeranian University of Technology, Szczecin, Poland, Head of the PhD Studies in the Electrical Engineering Faculty and the Head of the Department of Electrical and Computer Engineering. His research interests include electromagnetic field theory, analytical and numerical techniques and non-destructive testing of materials.

Adam Ryszard Zywica received MSc degree in Electrical Engineering from the West Pomeranian University of Technology, Szczecin, Poland, in 2014. Since October 2014, he is a PhD student at the Department of Electrical and Computer Engineering, Faculty of Electrical Engineering, West Pomeranian University of Technology, Szczecin, Poland. His research interests include non-destructive testing of materials.