Optimization of multi-angle Magneto-Acousto-Electrical Tomography (MAET) based on a numerical method

1.   Introduction

Dielectric characteristics of a tissue may reflect its physiological and pathological conditions. For instance, some cancerous tissues have been shown to present changes in electrical conductivity from the surrounding normal tissues ^[1,2,3], which potentially provides clinically relevant information for early diagnosis. Electrical Impedance Tomography (EIT) is a novel technique that measures conductivity with the advantages of high-speed, safety and low cost. However, its clinical application is as yet limited by poor spatial resolution due to the mathematically ill-posed problems ^[4,5]. Consequently, several novel methods have been proposed, such as Magnetic Resonance Electrical Impedance Tomography (MREIT) ^[6,7], Magneto-Acoustic Tomography with Magnetic Induction (MAT-MI) ^[8,9] and Magneto-Acousto-Electrical Tomography (MAET). MREIT is based, as its names suggests, on Magnetic Resonance Imaging (MRI). First, MRI images are acquired. Then, a weak electrical current is injected into a sample using electrodes. By using the J-substitution algorithm or the harmonic B_z algorithm, the image of local electrical conductivity can be mapped. MAT-MI was proposed by Bin He et al. in 2005 ^[10]. The electrical current in the test subject is created by magnetic field induction. Therefore, this method eliminates problems of discomfort and risks associated with current injection.

MAET is also known as Hall Effect Imaging (HEI) or the Lorentz Force Electrical Impedance Tomography (LFEIT) ^{[11,12,13,14,15,16]}. The method employs an ultrasound transducer to transmit a mechanical wave into the sample placed in a static magnetic field. As both positive and negative ions in the sample are displaced forward or backward along the wave propagation direction in the presence of the magnetic field, a distribution of Lorentz force is generated, and then an electrical field is created. By using one or more pairs of electrodes ^[12,17], placed around the sample, one or more voltage signals can be detected and analyzed. MAET was first demonstrated experimentally by Wen et al. ^[11]. In their study, 2D Hall effect images in an experimental setup were obtained. Then Montalibet ^[12,18] derived the measurement formula for MAET and used Wiener inverse filter to extract the conductivity parameter. Haider et al. ^[13] introduced the reciprocity theorem in the reconstruction algorithm. A high-spatial-resolution image of current density was obtained.

Grasland-Mongrain ^[14] further developed the model presented by Montalibet et al., and demonstrated that the detected voltage was proportional to the convolution of the electrical conductivity gradient with the acoustic pressure. Finally, reconstructed images of a gelatin phantom and a beef sample were presented. The relative amplitude of the front interface signal was different from the back interface due to the magnetic field inhomogeneity. Later, the concept of Acousto-Electrical speckle pattern was proposed by Mongrain et al. ^[15]. The theoretical similarity between the measured signal in LFEIT and ultrasound imaging modalities was assessed. Similar speckle patterns were observed, which was allowed to consider the use of ultrasound speckle-based image on LFEIT and to study electrical inhomogeneity structures.

Originated in radar applications ^[19], linear frequency-modulated pulse was applied to MAET by Sun et al. ^[16] in order to reduce peak stimulation power to the ultrasound transducer. B-scan images of various shaped agar phantoms were reconstructed finally. Dai et al. ^[20,21] studied multi-focus MAET method and obtained the conductivity parameter imaging of a rectangle phantom with a rectangle hole. Based on reciprocity theorem, Zengin ^[22] investigated the use of linear phase array (LPA) for MAET with magnetic measurements. The characteristics of the imaging system were analyzed by singular value decomposition (SVD) of the sensitivity matrix. The truncated SVD algorithm was implemented to reconstruct images with different signal noise ratio (SNR). Later on, the conductivity distribution in biological bodies was imaged by Gözü et al. ^[23] using LPA with electrodes measurements.

Most of the previous studies carried out experiments with linear scanning in which the conductivity boundaries in the sample were perpendicular to the ultrasound beam. This experimental setup may simplify the following image reconstruction because peak MAET signal occurs at the interface of two media with different conductivities when the interface is perpendicular to the ultrasound propagation path. However, for practical conditions with irregular interface, the MAET signal captured by the electrodes may decay and image reconstruction becomes challenging.

Recently, some researchers had tried to address this problem by obtaining more information from different scanning directions. Kunyansky et al. ^[17] proposed a rotational MAET method, with a rotating object of interest and two fixed pairs of electrodes. This method could reconstruct the boundaries of the regions of constant conductivity uniformly well, regardless of their orientations. However, the method was time-consuming by using the combination of 40 linear scan steps and 200 rotation angles.

Therefore, there is a need to improve the current method on both the accuracy of reconstructed image as well as reducing the total data acquisition time. In this study, a finite element model based on partial differential equation (PDE) of a multi-angle MAET experimental setup was constructed. In order to improve image resolution, as well as to lower computation cost, a higher ultrasound frequency was used. Consequently, the finite element method of Zengin ^[22] was substituted with a pressure acoustics time explicit interface method based on Range-Kutta discontinuous Galerkin (DG) method ^[24] in this study. To show the relation between the MAE signal and the conductivity, the corresponding measurement formula was demonstrated based on the study of Wen et al. Then, a model including two elliptically-shaped tumors was used. Next, a multi-angle MAET scanning method was proposed. Later on, multiphysics simulations using COMSOL Multiphysics software were implemented. The conductivity parameter images of complicated tissues were reconstructed by using multi-angle B-scan by image rotation. Finally, we investigated the relation between the number of angles and image quality and analyzed the signal-to-noise level based on the optimal number of angles. Compared with the method proposed by Kunyansky using similar scanning scheme and radon transform, our method had several advantages. First, the restriction on the angle of rotation was more flexible. Second, a simpler reconstruction was used, because the reconstructed method in this study was mainly based on image rotation. Then, the effect of the number of angles which was critical for the scanning time and the quality of reconstruction image was investigated in our study quantitatively with a large amount of numerical data. Finally, different signal-to-noise levels were investigated.

2.   Theory of Magneto-Acousto-Electrical Tomography

The model equations for the finite element simulation are described here. The geometry of MAET is shown in Figure 1. A medium with the electrical and acoustic properties of breast fat tissue with electrical conductivity spatial distribution $\sigma(y, z)$ is placed in a static magnetic field B₀ that is in the +x direction. A focus ultrasound transducer produces an ultrasound beam with particle velocity $\boldsymbol{v}(z, t)$ inside the breast fat. In order to simplify derivation, the acoustic attenuation due to viscosity is ignored for the propagation in the media. The Lorentz current density due to ultrasonically induced Lorentz fields is given by:

The Ohmic currents $\boldsymbol{J}_{O}(y, z, t)$ are related to the electrical potential $u(y, z, t)$ in the media by the Ohm’s law:

As the movement of charges is divergence-free, $\nabla \cdot\left(\boldsymbol{J}_{L}+\boldsymbol{J}_{O}\right) = 0$, we can obtain:

For the time explicit method used, the governing wave equations are formulated as the first order system, terms of the linearized continuity equation and the linearized momentum equation, as

where p, $\rho$ and c represent the acoustic pressure, the mass density of medium and the speed of sound, respectively. The acoustic pressure is obtained by the boundary condition of normal velocity on the transducer surface. The boundary condition is defined as:

where $v_{n}(t)$ is the vibration velocity produced by the transducer. The vibration velocity $\boldsymbol{v}(z, t)$ can be obtained by integrating against time,

Substituting (1.7) into (1.3) yields

The boundary conditions of the breast fat tissue are given differently on the electrodes and the other surfaces ^[25]. On the surface ($\partial \Omega$) without electrode:

and on each of the electrodes, the total current flowing out of the electrode is zero, since the electrodes are connected to a high impedance receiver:

where n is the outward normal on the boundary. The equations above are used in COMSOL to solve the forward problem of MAET.

Following Wen ^[11], the measurement formula of the MAET is given by:

where $\alpha, B_{0}$ and ${\rho _0}$ are constants. $R_{E}$ is the impedance of detection circuit. We can see that the MAE voltage is proportional to the gradient of conductivity. If the velocity of ultrasound and the conductivity gradient are in the same direction, the MAE voltage is the largest. On the contrary, if the velocity of ultrasound and the conductivity gradient are nearly perpendicular to each other, the MAE voltage will be approximately equal to zero. Hence, for irregularly-shaped tumors, rotating the sample is required to obtain an image of the complete pattern.

3.   Numerical studies

The COMSOL Multiphysics software (5.4, COMSOL Inc., Sweden) was used for the numerical studies. In this study, the pressure acoustics time explicit interface based on Range-Kutta discontinuous Galerkin (DG) method was used. Different from the pressure acoustics interface based on FEM, the pressure shape function used in DG was quartic (shape function in pressure acoustics interface was quadratic), so the mesh size in DG was set half of the wavelength to achieve appropriate spatial resolution. Consequently, the higher ultrasound frequency could be simulated with less solving time. The AC/DC (for electromagnetic) module of COMSOL Multiphysics was used to solve the forward problem of electromagnetic numerically. A numerical 2D geometry was built by COMSOL software. The distributions of focus ultrasound, current density and electric potential were obtained, respectively.

3.1.   Model geometry

The model geometry including normal tissue (breast fat) and malignant tumor (breast tumor) is shown in Figure 2. The flux density of the static magnetic field B₀ is 0.77 T along x axis. The breast fat is modeled as a 30 mm × 20 mm rectangle and breast tumors are modeled as two identical, 8 mm × 4 mm ellipses tilted at 15 and 75° off the z-axis. The tissue properties for sound speed, density and conductivity are listed in Table 1.

A focus ultrasound transducer with the radius of curvature of 50 mm is simulated. After the Helmholtz equation is solved, the frequency domain ultrasound pressure field is shown in Figure 3. The pressure distributions of z direction and y direction of are plotted in Figure 4a, b, respectively. The focus center is about 50 mm. Because the axis resolution of MAET is related to the duration time of pulse (code ultrasound is excluded), the sine pulse comprised of two cycle is adopted here. The input signal applied to the transducer is given as V(t) = 200sin(2πf₀t)H(t). The window function H(t) models the frequency response of an ultrasound transducer. Besides, the convergence of the result can be ensured by the window function. The voltage profile is depicted in Figure 5. According to the governing Eq (1.8) and boundary conditions (1.9) and (1.10), the forward problem of MAET can be solved by using COMSOL Multiphysics software. Note that the order of MAE signal is in μV level, so it is very necessary to set the relative tolerance less than 10e-5.

3.2.   Scanning process

In order to obtain B-mode image, the transducer is scanned along the y-direction with a step of 1 mm. The length of the tissue is 29 mm. After each B-mode scanning, the sample is rotated clockwise from 0 to 180° with the step angle of 180°/N. The N is the number of angles and it increases uniformly from 2 to 18 to evaluate the optimal N.

When the angle of rotation is 0°, the electrodes need to be placed to detect the induced electric field which is normal to the ultrasound, i.e. on the left and right sides of the sample. When the angles of rotation of sample are 0 to 40° and 140 to 175°, the electrodes do not need to change the position, and only need to rotate with the sample. However, when the angles of rotation of sample are 45 to 135°, the electrode needs to be replaced and placed on the other two sides of the sample at this time. There are two reasons for doing this, firstly to avoid the electrodes obstructing the ultrasound beam, and secondly to maximize the MAET signal.

3.3.   Current density

In the presence of a static magnetic field, an acoustic wave propagating in the sample induces a current density. Some researchers call it the velocity current density, since the distribution of the current density is consistent with the distribution of particle velocity of the longitudinal wave. When the focus transducer is placed in the y = 0 mm, the velocity current density at different times is shown in Figure 6. At this position, the ultrasound beams do not pass directly through the tumor regions with high conductivity. The region through which ultrasonic waves propagate is mainly homogeneous conductivity, so MAET signal will be induced at the front interface (z = 40 mm) and the back interface (z = 60 mm). And the voltage signal generated in the uniform conductivity region is minimal, which is consistent with all experimental results in other literatures.

When the transducer is placed in the position (at y = −6 mm), the beam propagates through the tumor, producing a large current density at the tumor position shown in Figure 7. The potential collected at both ends of the electrodes will contain the signal of the tumor area.

3.4.   Electric potential

The MAE signal is presented in Figure 8a, b when the transducer moves to 0 and −6 mm, respectively. By B-mode scanning along y axis, a MAE signal will be obtained on the electrodes after the transducer is excited using narrow pulse. Only the voltage signal of the front and the back interfaces (A and B in Figure 8a) of the sample can be detected on the electrodes when the transducer focuses on the non-tumor area. When the transducer focuses on the tumor area, not only the signal of the front and the back interfaces of the sample can be collected on the electrodes, but also the voltage signal of the tumor perpendicular to the ultrasonic direction can be obtained (A, B, C and D in Figure 8b). Note that in this model, the waveform of the MAE signal is similar to the ultrasound waveform transducer.

In this study, the distance between transducer and sample is 40 mm, and the sound velocity of water is 1490 m/s, so the data acquisition time of the first MAE voltage is 26.8 μs (A peak in Figure 8a). The width of the sample is 20 mm, and the sound velocity of breast fat tissue is 1454 m/s, so the time between the first interface and the last interface is 13.7 μs (B peak in Figure 8a).

3.5.   Data processing

In order to improve the quality of the signal, a Wiener filter is applied to the MAE signal. Then, the envelope of the MAE signal is obtained using the Hilbert transform, shown in Figure 9a, b. The subsequent B-mode images are reconstructed by using MAE signal from different excitation locations. Finally, after the linear interpolation algorithm is executed, the B-scan image using the relative amplitude of conductivity is obtained.

The total data acquisition time is 55 μs, and the data are sampled with time interval of 0.1μs. Thus, there are 551 samples in total. Considering the balance between image quality and processing time, we choose to perform eight times linear interpolation that results in 4401 numbers of sample. In the y axis, after linear interpolation, there are 1261 data totally. The reconstructed B-mode image of single angle is shown in Figure 10. The conductivity gradient along the ultrasound direction shows the greater intensity which is consistent with the theory.

In order to combine B-mode images with different rotation angles, we use the image rotation method in image processing. The B-mode image obtained at certain angle needs to rotate anticlockwise with the corresponding angle, and the rotated image is superimposed on previous image in a particular order. The reconstructed images with different number of angles are shown in Figure 11. The error analysis corresponding to the number of angles is demonstrated in section 3.6.

3.6.   Error analysis on the number of angles

It is expectable that more angles can improve the quality and accuracy of the reconstructed image. However, a tradeoff must be considered among the reconstructed accuracy, imaging time and complexity. We hope to obtain high quality images in reasonable time. In this section, error in the mean square sense is analyzed between reconstructed errors and the number of angles.

In this paper, we only evaluate the accuracy of reconstructed image profile, so we need to take a threshold on the reconstructed image to acquire corresponding profile. The ideal binary image based on the original conductivity distribution and practical reconstructed image are set as $\mu_{\mathrm{I}}$ and $\mu_{\mathrm{P}}$. Thus, we define

and quantify the difference between $\mu_{\mathrm{I}}$ and $\mu_{\mathrm{P}}$ in terms of mean-square-error (MSE) defined as follows:

The ideal conductivity parameter image should be given according to the length of the transducer emission waveform. The duration of the waveform is about 1.0 μs. The speed of sound in the tumor is 1700 m/s, so the length of the waveform is 1.7 mm. The resolution in theory is 0.85 mm. The ideal binary image is shown in Figure 12. The threshold of the reconstructed images is 0.25.

4.   Results

4.1.   Image reconstruction

Figure 11 presents the reconstructed conductivity parameter images using different the number of angles. It can be seen that when the number of angles increase from 2 to 10, the image quality is improved continuously, and the boundaries of the two tumors become more visible. By comparison, when the number of angles further increases from 12 to 18, the improvement of the image quality becomes negligible from the point of the reconstructed images.

4.2.   Optimal number of angles

Figure 13a shows the error data of the reconstructed image profile with respect to the number of angles. The values of MSE corresponding to the number of angles are listed in Table 2. When the number of angles increases, MSE decreases simultaneously. After the number of angles increases to 12, MSE is close to 0. It means under the corresponding threshold of 0.25, almost all profiles of irregular tumors are reconstructed with 12 angles, so 12 angles are more optimal considering the imaging time. The MSE curve was fitted by an exponential function $y = 1.632 e^{-0.392 x}$, with R² of 0.968. The threshold with respect to the MSE of 12 angles is depicted in Figure 13b. When the threshold is under 0.3, the accuracy of reconstructed image with 12 angles is acceptable.

4.3.   L² norm of the error with respect to measurement noise

In order to quantify the stability of the reconstructed method, we test L² norm of the error in the presence of the measurement noise. The signal-to-noise ratio (SNR) of measurement is defined as:$\mathrm{SNR} = 20 \log _{10}(\mathrm{S} / \mathrm{N})$. According to SNR, different noise levels (6dB, 14dB, 20dB, 40dB, 60dB and 80dB) are added to the MAE voltage to reconstructed images. In the experiment, there is usually a bandpass filter at the measurement front end, the MAE signal is processed by a bandpass filter (1–3 MHz) after adding noise. The reconstructed images with the measurement noise of (6dB, 14dB, 20dB, 40dB, 60dB and 80dB) are presented in Figure 14. When the noise level is 6 dB, the boundary of the reconstructed profile is contaminated by background noise, so the resolution of the image deteriorates. When the noise level is below 14 dB, the quality of reconstruction image is acceptable.

Because the noise is added to the reconstructed images, the MSE method with threshold is not fit to analyze the error here. Here, L² norm is introduced. The relative error is defined as:

where $\boldsymbol{x}_{R}$ and $\boldsymbol{x}_{O}$ are the reconstructed image and object image, respectively. Figure 15 shows the L² norm of the error with respect to the measurement noise.

5.   Discussion and conclusion

The present study employed numerical method to investigate some key aspects of a modified MAET imaging process. Most of the previous studies on MAET used linear scan with single element ultrasound transducer. To overcome the shortcoming of the conventional scanning method, this study proposed a novel multi-angle scanning method and implemented it in numerical simulation. The reconstructed image using the conventional scanning was shown in Figure 10. The boundary perpendicular to the transducer can be obviously observed while other boundaries are less visible or even invisible. In order to find optimal number of angles (related to step angle), 1044 scanning data were used with step angle 5°. The MSE shows that 12 angles were optimal considering the scanning time. The reconstructed images in the presence of the measurement noise were also illustrated. It means when the MAE signal was not seriously contaminated by noise, the reconstructed images were acceptable with our reconstructed method. Noise effect is found to be reasonably small as reconstruction error is low at SNR down to 25dB.

In order to improve the lateral resolution, the focus ultrasound is used in this study. The ultrasound beam width is not uniform in the longitudinal direction, and it may reduce the image area to keep the performance. Moreover, we have proposed a multifocus method in our previous study ^[20]. This method acquires multiple images with different focus depth and combines images to obtain the final image. Therefore, when the imaging area is larger, it is feasible to keep both the longitudinal and the lateral resolutions by using the multifocus image method.

The method proposed offers some advantages over the conventional methods. As shown in Figure 11, when more scanning angles are used, the detectable boundary becomes more completed clearly. Compared with the method proposed by Kunyansky ^[17] using similar scanning scheme, our method has several advantages. First, the angular range of our method is 0–180°, which is less than 0–360° of that method. Furthermore, this reconstruction method is more flexible. Thus, it may be used in MAET excited by linear phase array, because the deflection angle of linear phase array is difficult to cover the whole range of 360° in equal step. Second, the image reconstruction mainly based on the image transformation is very simple. Finally, our study quantitatively investigated the effect of rotation step which was critical for the scanning time and the quality of reconstruction image.

The data obtained from the present study have provided a proof of feasibility of a modified imaging procedure and reconstruction algorithm as well as an optimized configuration, which greatly reduces duration of data acquisition as well as ultrasound exposure. At present, an experimental study of a practical implementation of this method is ongoing.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant Nos. 61427806, 91859122, 81871429), and Shenzhen Basic Science Research JCYJ20170818141950351.

Conflict of interest

The author declares there is no conflict of interests.