Willmore-type variational problem for foliated hypersurfaces

1.   Introduction

Electrical impedance tomography (EIT) is a medical functional imaging technique that images the internal conductivity distribution by measuring boundary voltage data ^[1]. Under different physiological and pathological conditions, such as respiration, heart beating, cancer, the human conductivity and dielectric constant will change. EIT technology has been considered as an effective solution for human functional imaging, as it can reconstruct the electrical conductivity distribution of organisms. Compared to other tomographic techniques, such as computed tomography (CT) and magnetic resonance imaging, this technique has the advantages of being non-invasive, low-cost and radiation-free. It is widely used in biomedical imaging, such as lung ventilation monitoring ^[2], cancer detection ^[3] and brain imaging ^[4,5,6,7]. However, the image reconstruction of EIT is an ill-posed, nonlinear inverse problem ^[8], which causes EIT to suffer low resolution of the reconstructed images ^[9]. Therefore, it is of great value to develop the EIT image reconstruction method to obtain reconstructed images with clear boundaries, fewer artifacts and higher resolution images for clinical application.

Many methods have been proposed for EIT reconstruction. These methods can be divided into algebraic reconstruction techniques (ARTs) ^[10] and artificial neural networks (ANNs) ^[11]. ARTs reconstruct the conductivity distribution based on iterative back-projection. The regularization and iteration methods, including Landweber iteration ^[12], Gauss-Newton (GN) iteration ^[13], Tikhonov regularization ^[14] and total-variation regularization ^[15], have been implemented to reduce the influence of ill-posedness ^[16]. Nguyen et al. ^[17] proposed a time-efficient algorithm with a self-weighted NOSER-prior method and obtained results with good accuracy and noise tolerances. Sun et al. ^[18] presented an improved regularization method by combining the prior information extracted from the patient with Tikhonov regularization to obtain a high-resolution image. Jin and Maass ^[19] introduced the regularization theory with ℓp sparsity constraints to solve the inverse problems; Margotti ^[20] presented a strategy that combined a gradient-like method with a Tikhonov regularization to find stable solutions to ill-posed problems. Much research on ARTs has been done, but they still have problems such as a slow solving speed, blurry boundaries and vulnerability to noise.

Differing from the ART method, ANNs reconstruct images from voltage measurements through the use of a nonlinear transform, which requires even less computation when dealing with nonlinear problems ^[21]. However, when using deep learning methods for EIT image reconstruction, the spatial resolution of the reconstructed image is always lower than the training size. Therefore, the image reconstructed via EIT cannot be 100% restored, and this error will never tend to zero ^[22].

ANNs based on deep learning techniques have attracted much attention in the recent past for reconstruction of EIT images ^[23]. Wang et al. ^[24] developed a network by fusing the hybrid particle swarm optimization algorithm and radial basis function layer to obtain a compelling image. Huuhtanen and Jung ^[25] proposed an ANN with multilayer perceptrons and a radial basis function layer. However, these ANNs without convolutional calculation occupy too many memory resources, reducing the computational efficiency.

Therefore, the ANNs based on multilayer convolutional neural networks (CNNs) have been introduced to EIT analysis. Tan et al. ^[26] proposed an improved LeNet and proved the potential of CNNs to improve image accuracy and training speed. Gao ^[27] offered a CNN with a convolutional denoising auto-encoder and reconstructed a robust and denoising image. Hamilton and Hauptmann ^[28] combined the U-Net with the D-bar method to reconstruct both circle and lung phantoms. Although the method provided refined image development relative to the traditional D-bar algorithm, the method is highly dependent on the result of the D-bar algorithm. Wei et al. ^[29] used the improved U-Net composed of induced contrast current as the network input. They significantly improved the speed, stability and quality of EIT imaging, especially in sharp corners and edges. de Hoop et al. ^[30] and Nganyu et al. ^[31] also applied the concept of physics to solve partial differential equations and related parameter identification problems; it is also a feasible idea to solve inverse problems. Therefore, deep neural networks may become an effective method for image reconstruction based on EIT technology.

Because the exact inverse transform may not exist a priori, especially in the presence of sensor non-idealities and noise, EIT image reconstruction is a challenging task. To solve this problem, deeper and more complex deep learning models should be proposed ^[32]. The image reconstruction ability of deep learning models is reinforced with the increase of network layers, as they have multiple hidden layers that enable them to learn abstractions based on inputs ^[33]. But, deeper networks may cause gradient vanishing. Many deep learning models, including ResNet ^[34] and DenseNet ^[35], have been proposed to solve this problem. DenseNet is a meaningful choice to solve the EIT image reconstruction problem because it uses shorter connections between each layer of the dense block to maximize the information flow and improve the fitting ability. However, using the conventional DenseNet network still cannot improve the EIT image reconstruction task.

Considering the problems described above, a deep learning method based on DenseNet with multi-scale convolution (MS-DenseNet) is proposed here to improve the accuracy of EIT image reconstruction. To extract the information on different scale features, the proposed network uses convolutional kernels of different sizes to perform three types of convolutional operations, forming three types of dense blocks instead of the original dense block structure. The three dense blocks were placed in parallel to extract different input features. Additionally, aiming to minimize the information loss of the pooling module during image reconstruction, a hybrid pooling structure was designed and is applied in connection blocks. To solve the problem that the parameters do not easily reach the optimal solution during the training of large networks, a learning rate setting has been reduced to two stages applied to enhance the adaptive moment estimation (Adam) algorithm fitting ability. The proposed network was trained and tested on two datasets, including a circle-phantom dataset and a lung-phantom dataset. The voltage of the phantom was fed into the MS-DenseNet, and the predicted conductivity distribution was obtained via nonlinear calculation; the result was used to verify the method's effectiveness. The root mean square error (RMSE), structural similarity (SSIM), mean absolute error (MAE) and image correlation coefficient (ICC) were adopted as the evaluation metrics. The results showed that the MS-DenseNet method has good performance in terms of generalization and robustness as compared with the conventional DenseNet and the GN method; this helps to obtain EIT reconstructed images with clear boundaries and high resolution.

The contributions of this paper can be elaborated as follows:

1) The proposed MS-DenseNet method involves obtaining the conductivity distribution from the voltage inputs directly for EIT image reconstruction. The network incorporatess three types of dense blocks, a mixed pooling structure and a modified learning rate setting to boost the performance of the network.

2) Parallel dense blocks based on various types of multi-scale convolutions are used to replace the traditional dense block. Information from data of different scales is obtained and enhances the fitting and generalization capabilities of the network.

3) A hybrid pooling structure, which includes average pooling and 2 × 2 convolutions, has been used to replace the maximum pooling or average pooling in the conventional DenseNet, thus improving the information flow of voltage measurements and effectively preventing information loss.

4) A learning rate setting that is reduced to two stages is proposed to replace the fixed learning rate setting; it can significantly improve the fitting capability and optimize the image reconstruction results in the last relatively few epochs.

2.   Materials and methods

2.1.   Principle of EIT

The EIT image reconstruction process can be divided into forward and inverse problems. Figure 1 shows the process of EIT image reconstruction.

2.1.1.   Forward problem

The solution to the forward problem aims to obtain the electrode voltage ${\boldsymbol{U}}$ on the boundary $\partial {\boldsymbol{\Omega }}$ :

where $\sigma ({\boldsymbol{\Omega }})$ is a given distribution of conductivity in the region ${\bf{\Omega }}$ and $F(\cdot):\sigma ({\boldsymbol{\Omega }}) \to U$ is a forward mapping from $\sigma ({\boldsymbol{\Omega }})$ to ${\boldsymbol{U}}$. Let $\varphi ({\boldsymbol{r}})$ be the conductivity at any point ${\boldsymbol{r}}$ in the region ${\bf{\Omega }}$; the complete electrode model is used as forward mapping $F(\cdot)$; then, the forward problem can be calculated as:

where ${e_s}$ is the s-th electrodes set at $\partial {\boldsymbol{\Omega }}$, $n$ is the number of electrodes, ${z_s}$ is the contact impedance of electrode ${e_s}$, ${U_s}$ is the electric potential of ${e_s}$ and $v$ is a direction vector, where the direction is perpendicular to the tangential line at ${\boldsymbol{r}} \in \partial {\boldsymbol{\Omega }}$. The selection and the numerical solution of ${\boldsymbol{r}}$ in Eqs (2)-(5) are usually completed by using the finite element method (FEM) ^[36].

2.1.2.   Inverse problem

The solution to the inverse problem aims to reconstruct the image of the inner conductivity distribution ${\boldsymbol{σ }}$ in the region ${\bf{\Omega }}$ from the given boundary measurements ${\boldsymbol{V}}$. The traditional inverse problem can be described as minimizing the objective ${\boldsymbol{\Phi }}$ by finding the suitable distribution of conductivity ${\boldsymbol{\hat \sigma }}$ :

where $m$ is the number of voltage measurements at the boundary of the object, $\hat U$ is the measurement calculated by ${\boldsymbol{\hat \sigma }}$ and $G(\cdot)$ is a penalty function to prevent overfitting.

The deep learning-based EIT image reconstruction can be described as follows:

where $\zeta (\cdot)$ is the nonlinear mapping between ${\boldsymbol{V}}$ and ${\boldsymbol{\hat \sigma }}$, and the error between ${\boldsymbol{\hat \sigma }}$ and ${\boldsymbol{σ }}$ is reduced by modifying the parameters in $\zeta (\cdot)$.

2.2.   MS-DenseNet EIT image reconstruction

2.2.1.   Structure of MS-DenseNet

The overall structure of the proposed MS-DenseNet is shown in Figure 2. Based on the DenseNet, the proposed MS-DenseNet network consists of an input layer, dense blocks, transitional convolution-pooling layers and an output layer. Successive batch normalization (BN), rectified linear activation function (ReLU) and convolution (Conv) calculations are performed in dense blocks. Each dense block contains four BN-ReLU-Conv structures with the same output, interconnected by dense connections. The output of each dense block is the splicing of the input of that block with the output of the four BN-ReLU-Conv structures. Each transitional convolution-pooling layer connects two adjacent dense blocks, which contain a 1 × 1 convolution and a pooling layer. The size of the convolutional and pooling kernels varies according to the size of the weight matrix and the number of channels. The red arrow in Figure 2 indicates the direction of data flow, and the green arrow denotes the pooling layer. A set of three-dimensional vectors near each arrow represents the size of the data height, weight and number of channels. The number of final output channels "n" depends on the number of FEM elements in the conductivity distribution. Table 1 shows the details of each dense block, where "k" is the number of output channels of each BN-ReLU-Conv structure in the dense block. The number of channels of output data is equal to the number of filters in the final output module, independent of the number of channels in the input data.

The image reconstruction process for EIT is ill-conditioned and nonlinear. Therefore, acquiring and retaining more input features is the key to EIT image reconstruction. The design principle of MS-DenseNet is that the whole network can obtain and retain as much input feature information as possible, so as to yield higher-quality reconstruction images.

Dense blocks preserve input features through a densely connected structure. The role of multi-scale convolution in neural networks is to capture information from data at different scales. The combination of these two blocks makes the network more sensitive to input data, allowing it to obtain and preserve diverse features. The hybrid pooling can preserve more input features and effectively prevent information loss. The concatenation layer combines the features of the parallel input of the previous layer. These structures make MS-DenseNet suitable for EIT image reconstruction problems. The flowchart of the proposed MS-DenseNet is shown in Figure 3.

2.2.2.   Improved dense block structure

As shown in Figure 4, three different inception structures make up three kinds of dense blocks.

In the Convolutions type1 block, the output of the ReLU function is fed into the 1 × 1, 3 × 3 and 5 × 5 convolutional layers in parallel, and the number of channels corresponding to the filter is generated; then, the Sum layer sums the elements of the corresponding channels and corresponding positions in the output of different branches to obtain the computational results of the multi-scale convolution. Multi-scale convolution can yield feature information for different scales, which helps to improve the accuracy of reconstructed images.

In the Convolutions type2 block, the output of the ReLU function is first connected to a 1 × 1 convolutional layer to reduce the dimensions. The output of this 1 × 1 convolutional layer is then fed in parallel to a multi-scale convolutional layer. The Concatenation layer stitches together the features extracted by the three convolutional operations. Finally, the data dimensions are increased by the 1 × 1 convolutional layer. The connection layer combines the features of the parallel input of the previous layer to obtain more global information.

In the Convolutions type3 block, the output of the ReLU function is fed to the serial 1 × 1, 3 × 3 and 5 × 5 convolutional layers. The outputs of all three convolutional layers with different convolutional sizes are used as the inputs of the Sum layer. The result of the Convolutions type3 block is then obtained by the Sum layer. Various features extracted by multi-scale convolution use dense connections to reduce information loss during feature extraction.

In addition, 5 × 5 convolutions are replaced with 3 × 3 convolutions in Dense Block 3 due to the same calculation effect for the 2 × 2 input size.

2.2.3.   Hybrid pooling structure

In the EIT image reconstruction problem, the dimensionality of the output data is often larger than that of the input data, and any input data is valuable for the output. However, traditional pooling loses a large amount of data with unimportant features. Meanwhile, the convolutional layer can choose a proportional relationship between different inputs, allowing more detailed information to be retained.

As shown in Figure 5, the hybrid pooling structure combining the average pooling and 2 × 2 convolutions are used to replace the traditional max pooling or the average pooling structure in the connection block. To prevent information loss in the process of pooling, the feature maps of the hybrid pooling output are concatenated with the corresponding output feature maps of a 2 × 2 convolution and an average pooling layer. The average pooling layer provides 25% of the channels, and the 2 × 2 convolution provides 75% of the channels. The hybrid pooling structure can effectively improve information flow, reduce information loss and improve the quality of image reconstruction.

2.3.   EIT datasets

Two FEM phantoms, including a circle phantom and lung phantom, were used to complete the training of MS-DenseNet and DenseNet.

2.3.1.   Circle phantom

As shown in Figure 6, a circular tank with a radius of 9.5 cm was used as the imaging region for the circle phantom. The background conductivity was set to 1, two symmetrically distributed circles were selected to constitute the simulated normal conductivity region (blue region), the radius of the normal conductivity region was set to vary in the range of 3 to 4 cm, the normal conductivity region centers {x, y} were randomly distributed near the coordinates {5.5, 0} or {−5.5, 0} with a distance of less than 2.06 cm from the coordinates, the conductivity region was set to vary in the range of 0.1 to 0.5 and the step size was 0.05. Furthermore, 0-2 abnormal regions (red region) were set in the circle. The conductivity of the abnormal region was set to be 1.2 to 2, and the step size was 0.2. The radius of the abnormal region was set to 1-2.5 cm. Additionally, 20, 000 sets of data were established for network training to avoid overfitting. At the same time, it should be noted that the two abnormal conductivity regions cannot overlap.

The circle phantom datasets were taken from the EIDORS project ^[37]. Figure 6 shows that the process of creating the circle phantom datasets can be separated into three steps. The steps for creating the datasets are as follows:

1) Create a circle phantom, where the radius of the circle is 9.5 cm and 16 electrodes are evenly set in the circle's boundary. Divide the circle into 1650 elements via FEM generation.

2) Change the conductivity of normal region and abnormal region elements. Randomly generate 20, 000 sets of parameter data matrices, where each set of data includes center coordinates of abnormal regions, radius of abnormal regions, the number of abnormal regions, center coordinates of the normal region, and radius of the normal region. Furthermore, if an element is in both the normal and abnormal regions, its conductivity should follow the conductivity of the abnormal region.

3) Inject a current of 0.95 mA via a pair of adjacent electrodes; set the other two adjacent electrodes set as the voltage-measuring electrodes. Change the measurement electrodes until all of the electrodes (except for the injecting electrodes) are measured. Then, change the injected electrodes in sequence until all electrodes are injected. For each circle phantom with different conductivity distributions, collect 208 simulated measurements of voltage data as input data by using the above method; collect the conductivity of 1650 elements as output labels.

2.3.2.   Lung phantom

In order to simulate the monitoring of the lungs when the information is sufficient, a 2D lung phantom was created with actual lung and thorax boundaries. As shown in Figure 7, the lung phantom consisted of two lung regions and a background thoracic region. The conductivities of the lung and the thorax were set to 0.2 and 0.4, respectively, and the entire thorax region had the same conductivity. Three different conductivity anomaly types were applied in the lung region: Type a changes the conductivity of the whole lung region in the left or right, Type b changes the conductivity of a circle in the lung region and Type c changes the conductivity of a square in the lung region. Furthermore, the conductivity variation of the thoracic region was designed to be unaffected by these three variations. Considering the conductivity variation of the phantom in monitoring, the lung region conductivity was set to vary from 0.1 to 0.3, and the step size was set to 0.03. The thoracic region conductivity was set to vary from 0.35 to 0.45, and the step size was set to 0.02. All of the biological tissue conductivity values were obtained from a previous study about biological conductivity ^[38].

Chest CT images of an adult male were extracted from the EIDORS project and selected as lung and thorax segment points. The process of creating the lung phantom can be separated into the following three steps:

1) Load the sequential coordinates of different regional boundaries, including the lung and thorax. In this step, three sequential coordinates should be loaded, including a left lung, a right lung and a thorax.

2) Transfer the sequential coordinates to the FEM calculation software COMSOL Multi-physics to generate FEM meshes.

3) Use the created FEM meshes to solve the forward problem based on EIDORS data by using the same method as that for the circle phantom.

Each sample in the lung phantom dataset included 208 analog measurements of voltage data as input data and 3774 conductivity elements as output labels.

Both circle and lung data were randomly separated into two parts, including test datasets and training datasets. The test datasets included 1000 sets of data for the circular phantom and 1000 sets of data for the lung phantom. The training datasets included 18, 000 sets of data for the circular phantom and 10, 000 sets of data for the lung phantom. To further verify the stability and robustness of the proposed network structure, a 20 dB random Gaussian noise signal was added to the data of the test datasets to form datasets with 20 dB noise.

2.4.   Network training settings

For the proposed MS-DenseNet and DenseNet training, the labels were set by applying normalized conductivity distribution, including 1650 for the circle phantom and 3774 for the lung phantom. The sigmoid and the mean squared error (MSE) were selected as the output activation function and loss function, respectively, which are suitable for multidimensional regression tasks.

The network training was run on a PC with an Intel Core i7-9700, 3.00 GHz CPU and 16 GB RAM. The operating system of the computer was a 64-bit Windows 10, and the structure of the network was implemented by using the open-source deep learning library MatConvNet with MATLAB2019b implementation.

For the training process, 600 epochs with 500 mini-batches were set. The initial learning rates for the circle and lung phantoms were 0.0005 and 0.00005, respectively. The momentum was 0.9, and the L2 weight decay was set to be 0.001 to avoid overfitting. In addition, the learning rate was attenuated to 0.00005 for the circle phantom and 0.000005 for the lung phantom in the last 10 epochs. The loss function was MSE. The corresponding activation function was the sigmoid function. Adam is suitable for solving large-scale data and parameter optimization problems ^[39]. Thus, the solver used the Adam algorithm to keep training stable.

The conventional DenseNet and the iterative GN method were selected for comparison. The conventional DenseNet had the same training parameters as the MS-DenseNet, but the fixed learning rates were 0.0005 and 0.00005, respectively. And, they also trained with the same training datasets. The GN method was provided by the EIDORS project, using NOSER regularization. The maximum iterations were 10, the background conductivity was set to 1 for the circle phantom and 0.275 for the lung phantom and the hyperparameters were set to 0.55.

2.5.   Evaluation metrics

In order to compare the performance of different methods, RMSE, SSIM, MAE and ICC metrics were adopted to evaluate the reconstruction results. The RMSE, SSIM, MAE and ICC metrics between the reconstructed image and the actual image are respectively defined as follows:

where ${\sigma }_{i}$ and ${\widehat{\sigma }}_{i}$ denote the i-th element conductivity results of the reconstructed and the actual images, respectively, and n is the number of phantom elements. Additionally, ${\mu }_{\sigma }$ and ${\mu }_{\widehat{\sigma }}$ denote the mean of $\sigma$ and $\widehat{\sigma }$, respectively, ${\lambda }_{\sigma }$ and ${\lambda }_{\widehat{\sigma }}$ denote the variance of $\sigma$ and $\widehat{\sigma }$, respectively, and ${\lambda }_{\sigma, \widehat{\sigma }}$ denotes the covariance between $\sigma$ and $\widehat{\sigma }$. The parameters ${c}_{1}$ and ${c}_{2}$ are constants to maintain the stability of SSIM. L is the range of the conductivity value, ${k}_{1} = 0.01$ and ${k}_{2} = 0.03$.

3.   Results

3.1.   Training result

Figure 8 shows the average MSE between the test and training data according to epoch number for MS-DenseNet and the conventional DenseNet. Before the last 10 epochs, the learning rate settings of the two networks were the same. Figure 8(a) shows the case of the circle phantom. The MSE of the two methods quickly dropped to less than 5 in the first 30 epochs. Then, the rate of decline of the conventional DenseNet slowed down significantly and tended to stagnate after 300 epochs. The rate of decline of the conventional DenseNet is not apparent. MS-DenseNet maintained a significant decline after 300 epochs. At the same learning rate, compared with the conventional DenseNet, MS-DenseNet exhibited a better attenuation rate of the MSE curve; and, the final MSE was smaller during training. In the last 10 epochs, the learning rate of MS-DenseNet became one-tenth of the original. This part shows the comparison of MS-DenseNet before and after modification of the learning rate. In the last 10 epochs, due to the sudden decrease of the learning rate, the MSE finally declined further. The attenuation of the learning rate means that the step size of the search for the minimum MSE decreased. This adjusts the value of MSE closer to the minimum. Figure 8(b) shows the case of the lung phantom. Both methods maintained a faster rate of decline in the first 100 epochs and then maintained a steady and slow decline in the last 500 epochs. The reduction of the two methods became more unstable, especially the conventional DenseNet. It shows that the method in this paper yielded better results in terms of both the decay rate and the convergence value of the loss MSE curve given the same number of training epochs. The training time was about 117 hours, and the size of the parameters was 31.7 Mb.

3.2.   Noiseless dataset reconstruction analysis

The reconstructed image metrics calculated by MS-DenseNet, the conventional DenseNet and the GN method with noiseless test datasets are shown in Table 2. Compared with the conventional DenseNet, the proposed MS-DenseNet achieved an increase of 0.8, 1.5% (circle phantom), 0.4 and 0.4% (lung phantom) for the ICC and SSIM, respectively, as well as a decrease of 24.2, 28.6% (circle phantom) and 17.6, 18.7% (lung phantom) for the MAE and RMSE, respectively. Compared with the traditional GN method, MS-DenseNet achievevd a visible increase of 20.7, 25.3% (circle phantom) and 39.6, 41.6% (lung phantom) for the ICC and SSIM, respectively. The RMSE and MAE exhibited a visible decrease of 75.6, 65.9% (circle phantom) and 85.8, 89.4% (lung phantom), respectively. To facilitate comparative analysis, eight cases were randomly selected from the test datasets for observation.

Figure 9 shows the reconstructed images of eight cases extracted from the circle phantom test datasets. Table 3 shows the metrics of eight circle phantom cases on noiseless datasets. The reconstructed image of the GN method had an obscure edge and wrong estimates of the size and conductivity in the normal and abnormal conductivity regions. The implication of NOSER regularization in the GN method limited the sharp variety of conductivity distribution, which made the image edge smooth. The boundaries of different conductivity regions were difficult to distinguish. Furthermore, the region near the electrode was more prone to false predictions. The images reconstructed by the proposed MS-DenseNet and the conventional DenseNet revealed stable results. A clear and sharp edge was created between the background and the normal conductivity region in Case 2. The changes in the position and the size of the normal conductivity region were relatively small, so it was easier to reconstruct. The conventional DenseNet was unable to produce accurate images even in small abnormal regions (Cases 5 and 6); and, it yielded the worse reconstruction image in the central region (Cases 3-8). This is because the image reconstruction process of EIT is ill-posed and nonlinear, and because abnormal regions that are small or close to the center of the field will be more difficult to detect and reconstruct. Still, MS-DenseNet tended to predict two kinds of abnormal regions with a higher probability than the other methods. In Case 6, neither DenseNet methods could accurately predict the conductivity value in the region of abnormal conductivity. Compared with the DenseNet, MS-DenseNet predicted the presence of abnormal regions and yielded conductivity predictions that were closer to the truth.

Figure 10 shows the reconstructed images of eight cases extracted from the lung phantom test datasets. Table 4 shows the metrics of eight lung phantom cases on noiseless datasets. The GN method was unable to reconstruct the lung edge, and high conductivity artifacts were produced in the thoracic region. The two DenseNet-based methods accurately restored the boundaries of the lung conductivity region in Case 2. In Cases 1, 3, 4, 7 and 8, the conductivity of the anomaly region had a high variation. Both DenseNet-based methods were able to predict a change in lung conductivity; but, MS-DenseNet yielded a more accurate and more evident edge, while the conventional DenseNet underestimated the region of conductivity change or incorrectly estimated the conductivity value, as can be observed for the left lung in Cases 3 and 7. In Cases 5 and 6, the conductivity of the anomaly region had a small size, which is more difficult to predict. The conventional DenseNet could not recognize a conductivity change in the small region. However, the two DenseNet-based methods incorrectly estimated the conductivity value at the center of the variation region, which can be observed in the prediction of the left lung in Case 3.

3.3.   Noise dataset reconstruction analysis

In order to further verify the stability and robustness of different methods, test datasets with a 20-dB random Gaussian noise signal were used in the experiment. Table 5 shows the average metric results for 2000 datasets with 20-dB Gaussian noise. Compared with the noiseless datasets, the conventional DenseNet had a decrease in ICC and SSIM of 2.2 and 2.3% for the circle phantom and 0.9 and 1.2% for the lung phantom, respectively, as well as an increase in MAE and RMSE of 51.1 and 38.9% for the circle phantom and 45.9 and 53.6% for the lung phantom, respectively. The decline of four metrics in MS-DenseNet achieved 2.1, 2.1, 67.0 and 39.7% for the circle phantom and 0.6, 0.6, 59.0 and 40.0% for the lung phantom.

Figure 11 shows the same reconstructed images as Figure 9, calculated with the 20-dB signal-to-noise ratio (SNR) noise datasets. Figure 12 shows the same reconstructed images as Figure 10, calculated with the 20-dB SNR noise datasets. Table 6 shows the metrics of eight circle phantom cases on the 20-dB noise datasets. Table 7 shows the metrics of eight lung phantom cases on the 20-dB noise datasets. The GN method has almost no influence because stability comes from the NOSER regularization. The metrics of MS-DenseNet were still better. The 20-dB SNR noise did not have a significant influence on the normal conductivity region or the lung edge. Still, the conductivity prediction was more likely to lead to false distribution in the abnormal regions.

3.4.   Reconstruction of the image based on experimental measurements

An EIT system was used for practical verification experiments, as shown in Figure 13. Different concentrations of saline formed the background (1.0 S/m) or chest region (0.4 S/m), and their conductivity values were measured by using a conductivity meter. The normal region was replaced by nylon rods (0.0 S/m), and the agar model replaced the abnormal (2.0 S/m) and lung (0.2 S/m) regions. The device was composed of a measurement region, an analog acquisition board, a digital acquisition board and the upper PC. The circle image region was made of acrylic, where the inner diameter was 19 cm. (The lung region was the same size as the phantom.) The analog acquisition board had an AD8039 device as the voltage-controlled current source to release the excitation current. The voltage signal was connected to AD8251 and AD8253 devices to complete the amplification. An OPA1612 device was used as the voltage follower buffer. The digital acquisition board used a programmable digital signal processing development board to collect digital signals and control electrode conversion. The agar model was made by heating 6% agar powder, sodium chloride, 3% hydroxyethyl cellulose and 1% formalin solution. The conductivity of the agar model was measured by using a conductivity meter before solidification.

The network trained on the simulation training set was directly applied to the experimental models without training on the experimental data. A successful transition to the experimental data reflects the robustness of the proposed MS-DenseNet method.

Figure 14 presents the experimental reconstructed images for five cases. Table 8 shows the metrics for five cases under the conditions of experimental voltage measurement. The conventional DenseNet incorrectly judged the normal region in Cases 1 and 4. In addition, more artifacts were generated in the center of the image region in Cases 2 and 3. The results of the MS-DenseNet reconstruction were more stable than those of the other methods, especially in terms of the judgment of abnormal region.

4.   Discussion

The results of image reconstruction, shown in Figures 9-12, demonstrate the advantages of the MS-DenseNet method in terms of image accuracy. This section will discuss how certain operations or structures in the proposed method contribute to the performance of EIT image reconstruction. Each ablation experiment in the discussion used the same training datasets for experiments.

4.1.   Multi-scale convolution and different dense blocks

Multi-scale convolution can be traced back to GoogLeNet ^[41], where the inception block demonstrated the advantage of operating with multi-scale convolutional kernels instead of single-sized kernels. The role of multi-scale convolution in neural networks is to capture information from data at different scales, and obtaining more input features is essential for EIT image reconstruction ^[41]. For the same feature map, the calculation expectations of the single-scale and multi-scale convolution operation methods were the same. Still, the operation of the multi-scale convolutional kernel was able to process the feature information of different scales simultaneously, and the simple zero-padding operation was able to control the output data of different convolutional kernels to the same size. Then, the output data were spliced together so that the next layer could extract the feature information of different scales simultaneously, as well as improve the utilization of the network's computing resources. The role of parallel dense blocks is similar to that of parallel multi-scale convolution, which is to extract multiple kinds of features and fuse them.

In order to verify whether the new dense blocks can improve the network's performance, ablation experiments should be performed. Without changing any other structures, three different types of dense blocks were respectively replaced with conventional dense blocks to compare their MSE.

Table 9 shows the impact of different dense block types on the noiseless circle phantom final output. The results prove that all three dense block types can improve the network's fitting performance and reduce the MSE in the final output. In contrast to Dense Block type1 and Dense Block type2, Dense Block type3 applied 5 × 5 convolution after 3 × 3 convolution, which is equivalent to obtaining a larger-sized convolutional kernel ^[42]. Dense Block type3 reflects a better effect in this experiment, which is due to the stronger spatial correlation of the arranged data in the voltage matrix, and it provides a feasible theoretical basis for the future design of a specific convolutional kernel determined by the Jacobian matrix.

4.2.   Hybrid pooling

In the conventional DenseNet, the front layer uses maximum pooling to extract feature information, while the average pooling layer is used for subsequent connected blocks to preserve important high-dimensional details ^[35]. This is suitable for image classification and other problems. But, the role of the pooling layer in image classification and image reconstruction is different. In the case of image classification, the input image contains information that is not useful for determining the final label, such as background information. For EIT image reconstruction, any input voltage measurement will affect the final result, and highly correlated measurements will have a significant impact on elements with small changes. This means that input characteristics should be saved to the next layer as much as possible. Max pooling was inferior to average pooling, and average pooling was inferior to 2 × 2 convolution in terms of ability to preserve features ^[43].

As shown in Table 10, different pooling settings were implemented in the conventional DenseNet to search for better metrics of circle datasets without any noise. Max pooling yielded the worst effect among the three methods. This is because max pooling ignores most of the information and only keeps the maximum value in a 2 × 2 matrix. When average pooling was applied, the MSE decreased. The average pooling was affected by all pixels in the 2 × 2 matrix, thus retaining all valid information in the next layer. The average pooling compressed all 2 × 2 matrices in one dimension with the same weight, resulting in indistinguishable data compression for the network. The 2 × 2 convolution operation effectively adjusted the importance of data in different dimensions while keeping all data valid. Therefore, a small amount of 2 × 2 convolution can be added to the pooling layer channel to realize a substantial improvement. However, considering that the introduction of 2 × 2 convolution increased the storage of the network and the effect of increasing the number of 2 × 2 convolutions after 75% of the channel is not apparent, the final pooling structure consisted of 75% 2 × 2 convolutions and 25% average pooling layers.

4.3.   Learning rate setting

The largest initial learning rate was selected under the premise of ensuring the stable decline of the network, because the network uses Adam as the solver, which prevents severe sway in the later epochs of the network training compared to the stochastic gradient descent method. However, an appropriate learning rate setting is more conducive to the learning of neural networks. Four different learning rate settings were implemented in the conventional DenseNet with the noiseless circle phantom to find a suitable learning rate setting. The Type (a) setting reduced the learning rate to 0.1 times the initial learning rate after 300 epochs. Figure 15(a) shows that, after the learning rate was reduced, the MSE exhibited a significant decrease in the beginning; but, the reduction of the MSE then began to slow down, and the final result was not as good as before the modification. Similarly, the Type (b) setting reduced the learning rate to 0.1, 0.01, 0.001, 0.0001 and 0.00001 times the initial learning rate after 100, 200, 300, 400 and 500 epochs, respectively. As shown in Figure 15(b), the same significant decrease in the learning rate occurred in the beginning, but the final MSE values were still not good. In Figure 15(c), an impulse learning setting was used; it maintained the initial learning rate in most epochs, it was but 0.1 times the initial learning rate for the last five epochs of every 100 epochs. The final MSE was almost the same as that associated with the classic fixed learning rate, but every pulse was able to trap the MSE into a local minimum. In Figure 15(d), the implemented learning rate is as follows:

where s is the initial learning rate, $i$ is the number of epochs of the learning process and $s(i)$ is the learning rate of i epochs. The results show that a slowly decreasing learning rate will reduce the performance of the DenseNet on EIT image data. A large number of locally optimal solutions were distributed around each epoch of training, and a small learning rate caused the training to fall into a local optimum. However, the network can be trained to a local optimum after nearing the global optimum. Finally, the learning rate setting was determined as a fixed initial learning rate in the first 590 epochs and 0.1 times the initial learning rate in the last 10 epochs.

4.4.   Compared with other methods

To further demonstrate the performance of the proposed MS-DenseNet in terms of EIT image reconstruction, a comparative analysis was carried out with other methods that emerged in recent years. Table 11 shows a comparison of the evaluation metrics for EIT image reconstruction for different methods. The metrics in Table 11 are those based on noiseless datasets.

Although the evaluation indexes of different datasets cannot be directly compared, the datasets from the selected literature are similar to those of this study, which has certain reference value. These characteristics will not change when other similar datasets are used to train and test the network. The comparison results also verified the good performance of the proposed network in terms of fitting and generalization ability.

4.5.   Limitations of this study and future direction

MS-DenseNet was able to determine the exact locations and sizes of regions with different conductivities more accurately than traditional iterative methods. MS-DenseNet also has significant advantages over the conventional DenseNet, specifically in terms of the MAE and RMSE. But, in the case of the ICC and SSIM, the benefit was not as huge as imagined. The reason is that both the conventional DenseNet and MS-DenseNet were close to the maximum, which may be related to the complexity of the phantom design. In addition, when the size of the abnormal area was small and close to the middle area, the effect of the reconstructed image was worse than expected. In the future, the advantages of other traditional algorithms can be added to the neural network to further improve the robustness of the algorithm and use the dataset closer to the actual conductivity distribution of human tissue to enhance the network function.

With sufficient prior information, the method can image the conductivity distribution of human lungs and blood vessels in limbs. In the future, efforts on studying the problem of reduced noise immunity will be made, and the specific dataset will be more similar to the actual conductivity distribution of the human tissue to enhance the network function.

5.   Conclusions

This paper proposed an EIT image reconstruction method based on DenseNet with multi-scale convolution. Regarding the proposed network, it has three different multi-scale convolutions that form different types of dense blocks to extract feature information in parallel, which improved the generalization ability of the network. The hybrid pooling module, including a 2 × 2 convolutional layer and average pooling layer, was used to improve the information flow of the voltage data and reduce the loss of information during the pooling process. For the network training, a global learning rate setting has been proposed, which was demonstrated to optimize the fitting effect of the network. The simulation dataset and measured data were tested and analyzed with RMSE, SSIM, MAE and ICC as evaluation indicators. The results showed that the method proposed in this paper is competitive with other methods and can obtain reconstructed images with clear boundaries, fewer artifacts, higher resolution, high robustness and anti-noise ability.

Acknowledgments

This research was funded in part by the National Natural Science Foundation of China under Grant Nos. U22A20221, 61836011 and 71790614. The authors are also grateful for funding from the Fundamental Research Funds for the Central Universities 2020GFZD008, the 111 Project (B16009), Natural Science Foundation of Liaoning Province (2021-MS093, 2022-MS-119, 2021-BS-054) and the Basic Scientific Research Project of the Education Department of Liaoning Province in 2021(LJKZ0014).

Conflict of interest

The authors declare that there is no conflict of interest.