Three-Dimensional Holographic Electromagnetic Imaging for Accessing Brain Stroke

Abstract

The authors recently developed a two-dimensional (2D) holographic electromagnetic induction imaging (HEI) for biomedical imaging applications. However, this method was unable to detect small inclusions accurately. For example, only one of two inclusions can be detected in the reconstructed image if the two inclusions were located at the same XY plane but in different Z-directions. This paper provides a theoretical framework of three-dimensional (3D) HEI to accurately and effectively detect inclusions embedded in a biological object. A numerical system, including a realistic head phantom, a 16-element excitation sensor array, a 16-element receiving sensor array, and image processing model has been developed to evaluate the effectiveness of the proposed method for detecting small stroke. The achieved 3D HEI images have been compared with 2D HEI images. Simulation results show that the 3D HEI method can accurately and effectively identify small inclusions even when two inclusions are located at the same XY plane but in different Z-directions. This preliminary study shows that the proposed method has the potential to develop a useful imaging tool for the diagnosis of neurological diseases and injuries in the future.

1. Introduction

Medical imaging plays an essential role in the diagnosis of malignant tumors. Early diagnosis of the malignant tumor could significantly improve the treatment outcome and prognosis [1]. X-ray, ultrasound, computed tomography (CT), magnetic resonance imaging (MRI), and positron emission tomography (PET) are the commonly used biomedical imaging tools in hospitals. However, these techniques have some limitations. X-ray produces harmful radiation to the human body [2], and it is difficult for early abnormal tissue detection due to the relatively small contrast between the healthy tissue and the abnormal tissue at X-ray frequencies. PET is an excellent choice for imaging soft tissues, but it suffers from poor image resolution and relatively high-cost [3]. CT and MRI are well-established methods for identifying structural alterations of the biological tissue, such as brain tissue. However, they are unsuitable for continuously monitoring brain disease because of the relatively high-cost and time-consuming [4,5]. CT also produces harmful ionizing radiations to the human body. It is urgently needed to develop a new imaging method for biomedical imaging and diagnostic applications especially for early diagnosis of critical diseases, such as lung cancer, brain diseases.

Electromagnetic imaging (EMI) has been proposed as a potential technique to overcome the limitations of existing biomedical imaging tools, which has received increasing attention for applications in biomedical imaging and diagnostics. In the past three decades, several research groups have investigated EMI for imaging the electromagnetic properties (EPs, conductivity, permittivity, and permeability) of biological tissues using low-frequency (<10MHz) electromagnetic (EM) spectrum [6,7,8,9]. EMI has been applied to detect several diseases, including gastric emptying [10], lung ventilation [11], breast cancer [12], prostate cancer [13], pulmonary perfusion [14], and acute cerebral stroke [15]. EMI has been extensively studied for monitoring neurological diseases and injuries, including hemorrhagic stroke [16], ischemic stroke [17], and brain edema [18]. The applicability of EMI for monitoring brain diseases relies on the significant change in the EPs of brain tissues when a stroke occurs.

Electrical impedance tomography (EIT) is a noninvasive imaging tool which employs an array of electrodes on the surface of the biological object and maps the internal conductivity distribution changes of the object based on physiological and pathological activities [19,20]. However, the EIT technique requires expensive computational cost, may cause skin affections because of the use of the skin touch electrodes, and low spatial resolution. In 1993, Al-Zeibak et al. first reported the magnetic induction tomography (MIT) for biomedical applications [21]. An MIT system usually employs an excitation coil to induce eddy currents in the biological tissue and uses a detection coil to measure scattering field from the object to reconstruct a 2D image by using filtered back-projection. Compared to EIT, MIT is a contactless method which does not fill galvanic coupling in the space between the MIT device and the object. Therefore, the ill-defined electrode-skin interface is avoided. MIT is a more sensitive technique as it can monitor all three EPs parameters, and it is particularly attractive for identifying brain diseases (e.g., brain edema) because of the induced magnetic field via coils can penetrate through the skull much more easily than the injected currents via electrodes in the EIT system [22,23,24,25,26,27,28]. However, the MIT technique also has some drawbacks include limited image resolution, produce artifacts and the leak of clinical suitable measurement instruments. Recently, the authors developed a holographic electromagnetic induction (HEI) imaging for dielectric object imaging [29,30]. However, this method is unable to produce 3D images and unable to detect small lesions embedded in a 3D object accurately.

This paper presents the development of a 3D HEI method and investigates the feasibility and capability of the proposed method for detecting small inclusions embedded in a 3D dielectric object. A numerical system, including a realistic 3D head phantom and data processing model, has been developed to validate the effectiveness of the method for imaging infarcts that mimic a stroke embedded in the head phantom. This paper presents the theoretical model, together with simulation experiments and simulation results.

The present work is organized as follows: Section 2 presents the principle of 3D HEI imaging concept. Section 3 describes a numerical system for validation of the method. Section 4 presents simulation results. Finally, Section 5 concludes this paper.

2. Basic Principle

2.1. Concept of the 3D HEI System

Figure 1 shows a conceptual scheme of the 3D HEI system. The system consists of a cylindrical tank, a biological object, a radio frequency (RF) generator to produce EM signals, 16 excitation sensors to induce a magnetic field into the object, 16 receiving sensors to measure the scattering field from the object, a host computer with HEI program. All excitation sensors are equally arranged in a circle plane, and all receiving sensors are equally arranged in a circle plane. The excitation sensor array is located at the bottom of the tank, and the receiving sensor array plane is mounted on the wall of the tank and is designed moveable in the Z-direction. The biological object is located at the center of a cylinder (x = 0 mm, y = 0 mm, z = 0 mm). A matching medium is filled in the imaging chamber, which allows for an optimal sensor matching and ensures EM signals are propagating through the object and the scattering field can be more accurately recorded. For simplicity, a single frequency is selected as the working frequency according to previous studies [26,29].

We first place the receiving sensor array at one of the H vertical locations, the RF generator produces EM signals to excitation sensors, let one of the excitation sensors creates a magnetic field which leads to an electric field that drives an eddy current in the object, and every receiving sensor measures the backscattered signals from the object. This process is repeated until all excitation sensors have created a magnetic field. A 2D HEI image is obtained over the cross-section of the 3D object on the XY plane using 2D HEI method [29]. Then, we reposition the receiving array plane at a new vertical location and repeat the image data collection process until all vertical locations have been scanned. Vertical scans in the Z-direction can provide a stack of 2D images. A 3D object image can be reconstructed by combining a stack of 2D images.

2.2. Forward Model

This section presents the forward model that has been applied in the simulation experiments. The time-harmonic vector wave expression for the electric field can be solved by:

$$\nabla \times ({\mu }_r^{-1}\nabla \times \overrightarrow{E})-k_0^2{\epsilon }_r\overrightarrow{E}=-j\omega {\mu }_0\overrightarrow{J_s}$$

(1)

where $\nabla$ is the divergence operator, the relative permeability ${\mu }_r=\mu /{\mu }_0$, the complex relative permittivity ${\epsilon }_r=\frac{\epsilon }{{\epsilon }_0}+\frac{j\sigma }{{\epsilon }_0\omega }={\epsilon }^{\prime }+j{\epsilon }^{″}$, ${\mu }_0$ and ${\epsilon }_0$ are the permeability and permittivity of free space, respectively. $\sigma$ is the electrical conductivity, $\omega =2\pi f$, and $k_0=\omega \sqrt{{\mu }_0{\epsilon }_0}$, $\overrightarrow{E}$ denotes the electric filed, $\overrightarrow{E}={\overrightarrow{E}}_{inc}+{\overrightarrow{E}}_{scat}$ with ${\overrightarrow{E}}_{inc}$ and ${\overrightarrow{E}}_{scat}$ are the transmitted and received electric field in the sensor, $j=\sqrt{-1}$, and $\overrightarrow{J_s}$ is the excitation current density in the excitation coil.

The magnetic field strength is computed using the results of the electric field $\overrightarrow{E}$ in Equation (1):

$$\overrightarrow{H}=j{(\omega {\mu }_0{\mu }_r)}^{-1}·(\nabla \times \overrightarrow{E})$$

(2)

where $\overrightarrow{H}={\overrightarrow{H}}_{inc}+{\overrightarrow{H}}_{scat}$, ${\overrightarrow{H}}_{inc}$ and ${\overrightarrow{H}}_{scat}$ are the transmitted field and received magnetic field in the coil, respectively.

2.3. Backward Model

We assume a target point P embedded in a 3D object (see Figure 2). The sensor array is placed at one vertical location, the scattering signals from the object can be measured by any sensor located on the sensor array plane as [30]:

$${\overrightarrow{H}}_{scat}(\overrightarrow{r_0})=\frac{-j}{4\pi \omega {\mu }_0}\underset{V}{{\displaystyle \int }}[(\overrightarrow{J_m}·\nabla )\times \nabla +k_0^2\overrightarrow{J_m}+j\omega {\mu }_0\overrightarrow{J_s}\nabla ]G(\overrightarrow{r},\overrightarrow{r_0})dV$$

(3)

where $\overrightarrow{r_0}$ is the position vector from the origin to the point source, $\overrightarrow{r}$ is the position vector from the detector to point source, $G$ denotes Green’s function. The induced electric current density and magnetic current density can be computed by:

$$\begin{gathered} \overrightarrow{J_s}=j\omega {\epsilon }_0({\epsilon }_r-1)\overrightarrow{E}; \\ \overrightarrow{J_m}=j\omega {\mu }_0({\mu }_r-1)\overrightarrow{H}; \end{gathered}$$

(4)

Equation (3) can be simplified as:

$${\overrightarrow{H}}_{scat}(\overrightarrow{r_0})=\frac{k_0^2}{4\pi }\underset{V}{{\displaystyle \int }}[(a\overrightarrow{H}+b(\overrightarrow{H}·\widehat{\underset{\_}{r}})\widehat{\underset{\_}{r}})]G(\overrightarrow{r},\overrightarrow{r_0})dV$$

(5)

where $\widehat{\underset{\_}{r}}$ is the unit vector from the origin to the point source, $a={\mu }_r{\epsilon }_r-\frac{j({\mu }_r-1)}{k_{0}R}(1-\frac{j}{k_{0}R})$, $b=({\mu }_r-1)(-1+\frac{3j}{k_{0}R}+\frac{3}{{({\mathrm{k}}_{0}R)}^2})$, $R$ denotes the distance between the source and field point, $a$ and $b$ are proportional to $1/R^2$(i.e., $k_{0}R\ll 1$). Hence $k_0^{2}a\cong -({\mu }_r-1)/R^2$ and $k_0^{2}b\cong 3({\mu }_r-1)/R^2$.

The scattering signals from the region of interest with the negligible magnetic materials (${\mu }_r=1$) and with the non-negligible magnetic materials (${\mu }_r\ne 1$) can be rewritten as:

$$\{\begin{array}{}\mathrm{(6a)}& {\overrightarrow{H}}_{scat}(\overrightarrow{r_0})=\frac{-j\omega {\mu }_0}{4\pi }\underset{V}{{\displaystyle \int }}(\sigma +j\omega {\epsilon }_0{\epsilon }_r^{\prime }-{\sigma }_{bg})\overrightarrow{H}G(\overrightarrow{r},\overrightarrow{r_0})dV for ({\mu }_r=1)\hfill \mathrm{(6b)}& {\overrightarrow{H}}_{scat}(\overrightarrow{r_0})\cong \frac{1}{4\pi }\underset{V}{{\displaystyle \int }}\frac{{\mu }_r-1}{R^2}[-\overrightarrow{H}+3(\overrightarrow{H}·\widehat{\underset{\_}{r}})\widehat{\underset{\_}{r}})]G(\overrightarrow{r},\overrightarrow{r_0})dV for ({\mu }_r\ne 1)\hfill \end{array}$$

The tissue conductivities are relatively small at low-frequency RF spectrum, such as 10 MHz. In this case, the Born approximation can be applied in the above equation when performing predictive forward modeling of a given object and coil configuration [31]. The induced field inside the object can be modeled approximately as the incident field that exists at the same location without place the object in the imaging region.

2.4. Image Processing

Let us place the receiving sensor array at one of H vertical locations. The scattering signals from the target point P can be measured by each receiving sensor (see Figure 2). The visibility data of the scattering signals measured by a pair of receiving sensors can be calculated as [32]:

$$\overrightarrow{V_{vi}}(\overrightarrow{r_i},\overrightarrow{r_j})=\langle {\overrightarrow{H}}_{scat}(\overrightarrow{r_i})·{\overrightarrow{H}}_{scat}^{\ast }(\overrightarrow{r_j})\rangle$$

(7)

where ${\overrightarrow{H}}_{scat}(\overrightarrow{r_i})$ denotes the measured scattering field by receiving sensor located at $\overrightarrow{r_i}$, ${\overrightarrow{H}}_{scat}^{\ast }(\overrightarrow{r_j})$ is the conjugate complex of measured scattering magnetic field by receiving sensor located at $\overrightarrow{r_j}$, $<>$ denotes the time average. For a N-element receiving sensor array, the total visibility data is $\overrightarrow{V_{si}}={\displaystyle \sum }_i^N\overrightarrow{V_{vi}}(\overrightarrow{r_i},\overrightarrow{r_j})$, $i\ne j$, $N\ge 3$, at a selected vertical height (the distance between the object and the receiving sensor array). This visibility data contains phase and amplitude information.

Define the object density as [33]:

$$I(\overrightarrow{s})={\left(\frac{j\omega {\mu }_0}{4\pi }\right)}^2{\left|\sigma (s)+j\omega {\epsilon }_0{\epsilon }_r-{\sigma }_{bg}\right|}^2\overrightarrow{H}(\overrightarrow{s})·\overrightarrow{H^{\ast }}(\overrightarrow{s})$$

(8)

where ${\sigma }_{bg}$ is the conductivity of the medium.

The visibility formula over the object can be rewritten as:

$$\overrightarrow{V_{vi}}(\overrightarrow{r_i},\overrightarrow{r_j})={(\frac{j\omega {\mu }_0}{4\pi })}^2\underset{V}{{\displaystyle \iiint }}{\left|\sigma (s)+j\omega {\epsilon }_0{\epsilon }_r-{\sigma }_{bg}\right|}^2\overrightarrow{H}(\overrightarrow{s})·\overrightarrow{H^{\ast }}(\overrightarrow{s})\frac{e^{-j{k}_0(\overrightarrow{r_i}-\overrightarrow{r_j})·\underset{\_}{\widehat{s}}}}{s^2}dV$$

(9)

Equation (11) can be simplified by combining Equations (8) and (9).

$$\overrightarrow{V_{vi}}(\overrightarrow{D_{ij}})=\underset{V}{{\displaystyle \iiint }}I(\overrightarrow{s})\frac{e^{-j2\pi \overrightarrow{D_{ij}}·\underset{\_}{\widehat{s}}}}{s^2}dV$$

(10)

where $\overrightarrow{D_{ij}}$ is the baseline vector, $\overrightarrow{D_{ij}}=(\overrightarrow{r_j}-\overrightarrow{r_i})/{\mathsf{\lambda }}_{\mathrm{b}}$ with ${\mathsf{\lambda }}_{\mathrm{b}}$ denotes the wavelength of background medium, $\underset{\_}{\widehat{s}}=sin\theta cos\varphi \underset{\_}{\widehat{x}}+sin\theta sin\varphi \underset{\_}{\widehat{y}}+cos\theta \underset{\_}{\widehat{z}}$, $dV=s^{2}sin\theta d\theta d\varphi ds$. Change cartesian coordinates to spherical coordinates and define variables $p=sin\theta cos\varphi$, $q=sin\theta sin\varphi$, $n=cos\theta =\sqrt{1-p^2-q^2}$. Then, the element $dV$ can be represented as:

$$dV=s^{2}dpdqds/n$$

(11)

Substituting (11) into (10), visibility function becomes:

$$\overrightarrow{V_{vi}}(\overrightarrow{D_{ij}})=\underset{V}{{\displaystyle \iiint }}I(\overrightarrow{s})\frac{e^{-j2\pi \overrightarrow{D_{ij}}·\underset{\_}{\widehat{s}}}}{n}dpdqds$$

(12)

where $\overrightarrow{D_{ij}}$ can be represented by $u_{ij}=(\overrightarrow{x_j}-\overrightarrow{x_i})/{\mathsf{\lambda }}_{\mathrm{b}}$, $v_{ij}=(\overrightarrow{y_j}-\overrightarrow{y_i})/{\mathsf{\lambda }}_{\mathrm{b}}$, $w_{ij}=(\overrightarrow{z_j}-\overrightarrow{z_i})/{\mathsf{\lambda }}_{\mathrm{b}}$. If the coil array is positioned at a target vertical location, the visibility function can be expressed by the following equation:

$$\overrightarrow{V_{vi}}(u_{ij},v_{ij},w_{ij})=\underset{p}{{\displaystyle \int }}\underset{q}{{\displaystyle \int }}\underset{n}{{\displaystyle \int }}\frac{I(s,p,q)}{\sqrt{1-p^2-q^2}}{e}^{(-j2\pi {\Phi }_{ij})}dpdqds$$

(13)

where ${\mathsf{\Phi }}_{ij}=u_{ij}p+v_{ij}q+w_{ij}n$.

As all receiving sensors located on a planar plane ($w_{ij}=0$), the line integral can be obtained:

$$\tilde{I}(p,q)=\underset{s}{{\displaystyle \int }}\frac{I(s,p,q)}{\sqrt{1-p^2-q^2}}ds$$

(14)

The visibility function changes to:

$$\overrightarrow{V_{vi}}(u_{ij},v_{ij})={\displaystyle \iint }\tilde{I}(p,q)e^{-j2\pi (u_{ij}p+v_{ij}q)}dpdq$$

(15)

A 2D object image can be reconstructed by taking an inverse Fourier transform of Equation (15):

$$\tilde{I}(p,q)=\iint \overrightarrow{V_{vi}}(u_{ij},v_{ij})e^{j2\pi (u_{ij}p+v_{ij}q)}dudv$$

(16)

Equations (6), (7) and (16) can be applied to generate a 2D object image. The receiving sensor array plane is then placed at a new position in the Z-direction, and the 2D data collection process is repeated until the whole object has been scanned. Therefore, a 3D object image can be created by combining a stack of 2D HEI images.

As shown in Figure 2b, the depth information can be expressed by:

$$D_n=s_{n}cos({\theta }_n)$$

(17)

where ${\theta }_n$ is the radiation or detection angle of the sensor. Then:

$$dD=\frac{dz}{cos({\theta }_n)}=\frac{dz}{\sqrt{1-p^2-q^2}}$$

(18)

A 3D object image can be reconstructed by combining all 2D images that obtained when the sensor array repositioned at different vertical locations:

$$I(\overrightarrow{V_{vi}}=z_n,p,q)=\frac{d\tilde{I}(p,q)·\sqrt{1-p^2-q^2}}{dz}$$

(19)

The above equation can be represented as:

$$\frac{d\tilde{I}}{dz}=\frac{{\tilde{I}}_{z_n}+{\tilde{I}}_{z_{n-1}}}{z_n+z_{n-1}}$$

(20)

A 3D image can be reconstructed by combining Equations (19) and (20).

2.5. Metric

The peak signal-to-noise ratio (PSNR) can be applied to serve as objective criteria:

$$PSNR=10lo{g}_{10}(255/\sqrt{MSE})$$

(21)

where MSE (mean squared error) measures the quality of an estimator, $MSE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{(Y_i-\widehat{Y_i})}^2$, $Y$ is the vector of observed values of the variable being predicted, and $\widehat{Y}$ is a vector of $n$ predictions generated from a sample of $n$ data points on all variables. MSE is always non-negative, and MSE closer to zero demonstrates better performance.

The signal-to-noise ratio (SNR) measures the sensitivity.

$$SNR={\mu }_{sig}/{\sigma }_{bg}$$

(22)

where ${\mu }_{sig}$ and ${\sigma }_{bg}$ are the average signal value and the standard deviation of the background, respectively.

The structural similarity index (SSIM) measures the image similarity:

$$SSIM=\frac{(2{\mu }_x{\mu }_y+c_1)(2{\sigma }_{xy}+c_2)}{({\mu }_x^2+{\mu }_y^2+c_1)({\sigma }_x^2+{\sigma }_y^2+c_2)}$$

(23)

where ${\mu }_x$ and ${\mu }_y$ are the average of the original image $x$ and reconstructed image $y$, respectively. ${\sigma }_x^2$ and ${\sigma }_y^2$ are the variance of $x$ and $y$, ${\sigma }_{xy}$ is the covariance of $x$ and $y$. The two variables $c_1={(0.01L)}^2$ and $c_2={(0.03L)}^2$ to stabilize the division with the weak denominator, $L$ is the dynamic range of pixel-values. Where $-1\le SSIM\le 1$, and $SSIM=1$ can be obtained when $x=y$.

3. Numerical Experiments

A numerical system, including a realistic head model, a 16-element coil array, and imaging process model, was developed to evaluate the effectiveness, sensitivity, and accuracy of the proposed method for imaging infarcts that mimic a stroke embedded in the realistic head model. All simulations were conducted using MATLAB 2018a (The MathWorks, Inc. Natick, MA, USA) on a Dell Precision 5820 workstation which has an Intel Xeon W-2145 CPU with an internal clock frequency of 3.7 GHz and 256 GB of memory.

3.1. Simulation Setup

Figure 3 illustrates the top view of the numerical system setup. The system consisted of 16 excitation coils and 16 receiving coils mounted on the wall of a cylindrical tank (80 mm in radius, 100 mm in height). All excitation coils were equally arranged in a circle, which positioned under the cylindrical tank. All receiving coils were equally arranged in a circular array, and this array plane was designed moveable in the Z-direction. A 3D head phantom was placed at the center of the cylindrical tank (x = 0 mm, y = 0 mm, z = 0 mm). A high dielectric medium (${\epsilon }_r=88$, $\mathsf{\sigma }=0.48$ S/m, similar to distilled water) was filled in the imaging chamber (cylindrical tank), and a single frequency of $10$ MHz was selected as working frequency based on previous studies [26,29]. To collect 3D image data, the head model was placed at z = 0 mm, the sensor array place was moved from axial z = −50 mm to z = 50 mm in steps of 1 mm. For each vertical location, a 2D image data set was recorded to produce a 2D image. This data collection process was repeated until the entire head model has been scanned in Z-direction and a total number of 3,840,000 3D image data was collected. A 3D object image was reconstructed by combining a stack of 2D images.

3.2. RF Coil

For simplicity, the single circle loop coil (radius of 25 mm) was simulated as excitation sensor and receiver [29]. The magnetic field induced by each excitation coil can be computed by using Biot-Savart law:

$$\overrightarrow{H}=\frac{{\mu }_{0}I}{4\pi r_c{}^2}\oint dL=\frac{{\mu }_{0}I}{4\pi r_c{}^2}2\pi r_c=\frac{{\mu }_{0}I}{2r_c}$$

(24)

where $I$ is the current in excitation coil, $r_c$ denotes the radius of the loop coil, and $L$ denotes the length of the coil. The current in the coil was 1 A.

3.3. Head Model

An ellipsoidal shaped head model (long radius 60 mm, short radius 55 mm, thickness radius 40 mm) was developed to provide the validity of the proposed method. The head model consisted of a skin layer (3 mm thick), fat (5 mm thick), skull (7 mm thick), cerebro-spinal fluid (CSF, 3 mm thick), grey matter, white matter, and dura. Sphere-shaped inclusion was inserted into the head phantom to simulate infarcts that mimic a stroke. The brain lesion was assumed to consist of 100% blood. In this study, two inclusions (different sizes and positions) have been investigated. We used the published EPs values of the tissues to develop the head model [34].

Table 1. EPs of tissues at 10 MHz [34].

Structure	Thickness (mm)	σ (S/m)	Relative Permittivity
Matching medium	20	0.48	88
Skin	3	0.1973	361.66
Fat	5	0.029	13.767
Skull	7	0.0828	53.8
CSF	3	2	109
Grey matter	6	0.292	320
White matter	7	0.159	176
Dura	5	0.544	194.93
Blood	3/5	1.097	280

shows EPs of tissues at a frequency of 10 MHz. To save the computational time, the head model contained $80\times 80\times 50$ voxels which correspond to a region of $160\times 160\times 100$ mm³. The EPs of tissues have been assigned for each voxel to develop head phantom. The Gaussian function was applied to demonstrate the head model and brain lesion model.

4. Results

Several simulations have been performed to test the detectability of inclusions embedded in the multilayer head model using the proposed method. The first simulation performance was carried out with the conductivity set to zero and the relative permittivity set to unity (${\epsilon }_r=$ 1) for all tissues, thus simulating the host medium. Figure 4a,b display the real part (relative permittivity) and imagery part (conductivity) of the model under test, and its 3D reconstructed images are shown in Figure 4c,d.

The second simulation experiment was conducted with the conductivity and relative permittivity values set to the correct information (see Table 1), with no brain lesion present. Figure 5a,b show the real part (relative permittivity) and imagery part (conductivity) of the head model with no brain lesion and its 3D reconstructed images are demonstrated in Figure 5c,d. Results show that the internal structures of the head model can be clearly identified in the imagery part of the 3D reconstructed image. However, not all internal structures of the head model can be observed in the real part of the 3D reconstructed image. The computational time was 303.98 seconds for the first two simulation experiments.

The third experiment was carried out with the head model contains two lesions present in the imaging chamber. As shown in Figure 6a,b, the head model under test contains two lesions where brain lesion A (5 mm in radius, squared in black) located at (0 mm, 0 mm, 8 mm) and brain lesion B (5 mm in radius, squared in black) located at (0 mm, 0 mm, −8 mm). The real part and imaginary part of the 3D reconstructed images are shown in Figure 6c,d. Results show that two inclusions can be clearly identified in the imagery part of the 3D reconstructed image with the correct size, shape, and location information. However, no inclusion can be observed in the real part of the reconstructed image. The total computational time for this simulation experiment was 297.232 s.

The fourth experiment was carried out with the head model contains two different size lesions present in the imaging chamber. As shown in Figure 7a,b, the head model under test contains two lesions, brain lesion A (3 mm in radius, squared in black) located at (0 mm, 0 mm, 0 mm) and brain lesion B (5 mm in radius, squared in black) located at (0 mm, 0 mm, −8 mm). The real part and imaginary part of the 3D reconstructed images are shown in Figure 7c,d. Results show that two inclusions can be clearly identified in the imaginary part of the 3D reconstructed image with the correct size, shape, and location information. However, no inclusion can be observed in the real part of the reconstructed image. The total computational time for this simulation experiment was 297.232 s.

We have compared the proposed method with the developed 2D HEI method. To collect 2D image data, we have placed the sensor array plane and the head model at z = −50 mm and z = 0 mm, respectively. Figure 7e,f present the real part and imagery part of the 2D reconstructed images (XY plane) of the head model, respectively. However, only one of two inclusions can be observed in the imagery part of the 2D reconstrued image. The total computational time to produce the 2D images was 157 s.

The fifth simulation was conducted using the head model contains two inclusions (5 mm in radius, squared in black) located at (10 mm, 10 mm, 6 mm) and (0 mm, 0 mm, −6 mm). Figure 8a,b show the permittivity and conductivity values of the head model. Figure 8c,d display the real part and imaginary part of 3D reconstructed images, respectively. Figure 8e,f present the real part and imaginary part of the 2D reconstructed images, respectively. Results show that two inclusions can be observed in the imaginary part of the 3D reconstrued image with the correct size, shape, and location information. However, only one of two inclusions can be observed in the imagery part of the 2D reconstrued image.

Color bars in the original images plot the permittivity and conductivity values of the model under test, while color bars in the constructed images demonstrate the backscattered energy density distributions in the reconstructed images. Not all internal structures of the brain tissues can be identified in the reconstructed images especially the real part (permittivity) of the reconstructed images. The conductivity reconstructions (imaginary part) have higher quality than the permittivity reconstructions (real part), this finding agrees with the previously published research outcomes [35,36]. A total number of $(16\times 16\times 15=3840)$ 2D image data, and ($16\times 16\times 15\times 100=384,000)$ 3D image data were measured to produce a 2D image and a 3D image using the proposed method. Thus, more helpful information can be obtained in the reconstructed 3D images compared to the 2D images.

Table 2. Simulation errors.

Model	Head Contains Two Inclusions (5 mm in Radius) Located at (0 mm, 0 mm, 8 mm) and (0 mm, 0 mm, −8 mm)	Head Contains Two Inclusions (3 mm, 5 mm in Radius) Located at (0 mm, 0 mm, 0 mm) and (0 mm, 0 mm, −8 mm).	Head Contains Two Inclusions (5 mm in radius) Located at (10 mm, 10 mm, 6 mm) and (0 mm, 0 mm, −6 mm)
PSNR (dB)	49.9943	36.6764	49.9943
SNR (dB)	5.6384	2.0651	5.6382
MSE of real part	0.9657	0.9657	0.9657
MSE of Imagery part	0.4896	0.0876	0.0498
SSIM of real part	0.0455	0.4896	0.04986
SSIM of Imagery part	0.0455	0.0455	0.0455

shows simulation errors. It can be seen that with the proposed method it is much easier to identify larger size inclusions. Compared to previous published 2D image results [29], the proposed method could detect two inclusions which were located at the same XY plane but in different YZ planes.

5. Conclusions

This paper reported the theoretical framework of 3D HEI for internal structure imaging of small inclusions embedded in biological tissue. A numerical system, including a realistic 3D head model, a 16-element excitation sensor array, and a moveable 16-element receiving sensor array as well as image processing model, has been developed to test the detectability of lesion using the proposed method. Several simulation experiments have been conducted to evaluate the effectiveness, accuracy, and sensitivity of the proposed framework for detecting brain lesions embedded in the head model. A comparison study between the 2D HEI and 3D HEI has been conducted through various simulation experiments. The results showed that the 3D HEI could produce 3D brain images and detect inclusions with the correct size, shape, and location information. The two inclusions could be successfully observed in the reconstructed 3D images even when they were located at the same XY plane but in different YZ planes.

The HEI image quality highly depends on the visibility data that compares the scattering signals measured by any pair of receivers. For 16 excitation sensors and a 16-element receiver array that placed at 100 different vertical locations, a total number of $\left[16\times 16\times (16-1)\times 100=384,000\right]$ data can be recorded to produce a 3D image using the proposed method, while only a total number of ($16\times 16=256$) data can be measured to produce a 2D image using the conventional MIT method. Thus, the proposed method has the potential to produce a higher resolution image compared to the conventional EM techniques. Previous experimental studies [37] have shown that the 2D image quality also depends on the sensor array configuration and the sensor number. However, whether this claim is suitable for the practical 3D HEI system is one of the target future works.

Head is a very complex organ with functions and mechanisms have not been fully discovered. EMI technique exploits structural data of different levels, such as EPs of brain tissues, will help one to understand the internal structure of brain tissues, which offers the potential for detecting and monitoring brain disease. However, only ideal models (noise free) were involved in this study; this may affect the effectiveness of the proposed method. The image quality could be degraded due to the noise acquired from the practical measurement instrument. Furthermore, the proposed method requires further validation in more realistic experiment scenarios. From an experimental point-of-view, develop an EM device to collect complete EM data for biomedical applications remains both highly desirable and extremely challenging.

The obtained theoretical results showed that the proposed framework has the potential to become a useful tool for realizing complex situations and optimizing the whole system before the device is tested in practice. Future work should evaluate more realistic experimental data gathered from physical models, biological tissues, animals, and human subjects, to ensure that the proposed method is robust to more realistic scenarios. The development of a sensitive sensor for practical implementation of the proposed method is another target future research direction.