Deep learning for image reconstruction in thermoacoustic tomography

The neural network used for image reconstruction is constructed based on the encoder–decoder structure. The corresponding neural network architecture of the TAT-Net is illustrated in Fig. 1. In the TAT-Net, the encoder part is responsible for the image downsampling, which continuously extract the image features to obtain the advanced feature maps. Conversely, the decoder part implements continuous up-sampling, which maps the extracted features to the final reconstructed image. The input of the TAT-Net is the two-dimensional (2D) sinogram which is processed by arranging the sensor data collected in all directions (i.e., sinogram in a size of 270 × 180 represents that the number of channels is 180 and the data length of each channel is 270), and the output is the initial energy density distribution. The input image is fed into a 7 × 7 convolution layer with 2 × 3 stride to obtain the square feature map. The encoder part consists of 12 blocks, including 8 B1 blocks and 4 B2 blocks, and the decoder part also consists of 12 blocks, including 7 B1 blocks and 5 B3 blocks. The pooling layer and up-sampling operator are replaced by convolution layer and transposed convolution layer with a stride of 2 × 2 for higher representation capacity. To avoid the checkerboard artifacts,^[58] the size of the transposed convolution kernel is set to be twice the stride.

Zoom In Zoom Out Reset image size

Unlike previous neural networks for reconstruction task such as U-Net, the global residual learning and feature fusion between encoder and decoder are not recommended in TAT-Net. There exists a significant difference between the input sinogram and output reconstructed image, and learning the residual between them is more difficult, as using global residual in this task leads to a high probability of unconvergent training. The shallow feature fusion will result in shallow features being output directly into the final image without sufficient convolution for feature translation, which introduces some sinogram-like artifacts into the reconstructed image (see Subsection 3.2). In addition, we observe porphyritic BN artifacts^[59] in the experiments. Thus, we remove all batch normalization layers from residual blocks for stable training and consistent performance.

Since we remove the batch normalization layers, training a deep network becomes difficult. We adopt the residual module^[60] to solve the problem of gradient disappearance. To reduce training time and memory consumption, bottleneck structure is used in the entire TAT-Net. The B1 block is a standard bottleneck residual block. In the bottleneck residual structure CNN, traditional up-sampling and down-sampling operations take up a large part of the overall time and computation of the network, so we design the downsampling blocks B2 and up-sampling blocks B3 with bottleneck residual structure. The residual shortcut used in the down-sampling block is a 2 × 2 stride convolution and up-sampling is 2 × 2 stride transposed convolution.

The TAT-Net uses leaky rectified linear units (a = 0.01) as activation function to obtain a steady training. We initialize the weights in convolution layers and transpose the convolution layer into a Kaiming normal distribution (a = 0.01).^[61] The parameters of the network are optimized by minimizing the mean square error (MSE) loss function:

$$\begin{eqnarray}&&L={\Vert GT-f(X)\Vert }_{2}^{2},\end{eqnarray} \tag{ 1 }$$

where GT is the ground truth, X is the sinogram image, and f is TAT-Net operator.

We continuously reduce the MSE loss between the model output and GT by Adam algorithm^[62] to optimize the parameters of the model. The learning rate and batch size set to be 10⁻⁴ and 8, respectively. In the prediction process, we input the sinograms into the trained model to obtain the initial energy loss density images.

In the experiments, the network is implemented in PyTorch framework.^[63] and deployed in a computer platform with Intel(R) Xeon(R) CPU E5-2630 v4 2.2 GHz and NVIDIA TITAN Xp with 12G memory.

The deep learning-based method depends on a large quantity of data; however, the current experimental conditions of TAT do not have the capacity for data mass production. In addition, the data and images acquired through the experiments can be affected by various factors, such as system instabilities, difference among the reconstruction methods, and inhomogeneous distribution of the microwave field, etc. These issues will impede us from obtaining realistic target parameters, and thus, we seek for a reliable simulation method to construct the synthetic dataset that satisfies the mapping of the real initial energy density distribution into the sensor data.

Our data generation method includes two steps. The first step is to solve the forward solution of the thermoacoustic wave equation through the finite element method (FEM),^[64] and the thermoacoustic wave equation is defined as

$$\begin{eqnarray}&&{\nabla }^{2}p({\boldsymbol{r}},t)-\displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}p({\boldsymbol{r}},t)=-\displaystyle \frac{\beta \varphi ({\boldsymbol{r}})}{{C}_{{\rm{p}}}}\displaystyle \frac{\partial J(t)}{\partial t},\end{eqnarray} \tag{ 2 }$$

where c is the acoustic speed, β is the volume expansion coefficient, C_p is the specific heat, φ ( r ) denotes the energy deposition in the tissue, and J(t) represents the microwave pulse. Equation (2) is based on the assumption that microwave power dissipation in a region is so rapid that thermal diffusion can be neglected. J(t) = δ (t – t₀) is also assumed in our study.

In the FEM simulation, a total of 180 acoustic receivers are equally distributed along the surface of the 80-mm diameter circular background region. The finite-element mesh used for the forward calculations has 5371 nodes and 10560 elements. One to three circles of different sizes and energy densities are deployed as initial energy densities, and these targets are randomly placed at different positions in the region of interest (ROI) with 40 mm × 40 mm in size. We assume that the distribution of ultrasonic velocity and microwave energy are uniform in simulation. We can obtain sensor data S^c (N) from each acoustic receiver's channel, where N is the data length.

The second step is to preprocess the data. In order to reduce the computational cost of the neural network, we reduce the input size of the network by deleting redundant data. According to the size of ROI and the location number of the acoustic receivers, we can calculate a data range that contributes significantly to the reconstructed image:

$$\begin{eqnarray}&&{S}_{{\rm{new}}}^{c}={S}^{c}(i),\,i\in \left(\displaystyle \frac{{f}_{{\rm{s}}}({R}_{0}-a)}{\nu },\displaystyle \frac{{f}_{{\rm{s}}}({R}_{0}+a)}{\nu }\right),\end{eqnarray} \tag{ 3 }$$

where f_s is the sampling frequency, R₀ is the radius of the circular acoustic receiver array, a is the size of square ROI, and ν is the ultrasonic velocity. We use the following parameter configuration, f_s = 10 MHz, and ν = 1495 m/s. In addition, Gaussian noise and random bias are added into the sinogram to improve the ability to generalize the model and to reduce the domain gap between synthetic and realistic data. The processed sinogram has a signal-to-noise ratio (SNR) of 20 dB.

Finally, we use 270 × 180 pixel sinogram as the network input and 192 × 192 pixel energy density image as the ground truth. The dataset consists of 8192 training sets and 1024 validation sets. In order to verify the robustness of the model in terms of the number of imaging targets, an additional data containing 4 or 5 circular targets are generated.

To further demonstrate the feasibility of our proposed method, we conduct several phantom experiments with a TAT system shown in Fig. 2(a) to obtain the phantom experiment data. Plastic tubes containing physiological saline solution are used as the reconstruction target, and the photograph of the three tubes is shown in Fig. 2(b). In the TAT system, a microwave pulse is generated by a custom-designed 3-GHz microwave generator (bandwidth: 50 MHz, peak power: 70 kW, pulse duration: 750 ns), couples into a rectangular waveguide and is irradiated to a phantom by a horn antenna (114 mm × 144 mm). The averaged microwave power density at the phantom surface is less than 0.32 mW/mm², which is far below the safety standard (20 mW/cm² at 3 GHz). A single-element unfocused transducer (V323, Olympus, central frequency: 2.25 MHz, nominal bandwidth: 1.35 MHz) is used to acquire the TA signal over 360° at 180 positions around the sample rotated by the step motor. The digital signals are amplified by an amplifier and then recorded by a data acquisition card (PCI4732, Vidts Dynamic, sampling frequency: 50 MHz). The data are down sampled from 50 MHz to 10 MHz and the ROI is the same as that in the simulation setup. In the experiment, the scanning radius of the probe is 164 mm. In order to make the experimental data and simulation data have the same scanning radius property, a method of processing the signal is used to ensure the same detector positions as the simulation setup by setting the virtual detectors.^[65] As with the processing of synthetic data, we perform data refinement on the phantom data according to Eq. (3). Since the reflected signals created by non-target components in a TAT system factor are often present in these deleted data, the data refinement has the additional benefit of avoiding misinterpretation of the network due to the reflected signals. So far, we have obtained 270 × 180 pixel thermoacoustic experimental sinogram as the test dataset.

Zoom In Zoom Out Reset image size

We evaluate the performance of the proposed method by comparing the reconstruction results of DAS, U-Net, and TAT-Net. Figure 3 shows the results from the synthetic data with 1 to 3 circular targets. We can find that the TAT-Net can accurately reconstruct the position and shape of each target. Comparing with DAS, the image from TAT-Net has more well-defined boundaries and purer background. The U-Net provides more satisfactory results than DAS but produces additional artifacts in the 3-circle case. To visually show the reconstruction of TAT-Net, in Fig. 4 are plotted the one-dimensional (1D) line intensity profiles passing through the red dashed lines shown in Fig. 3 for the images reconstructed by different methods. It is obvious that the line intensity profile from our proposed method is closer to the ground truth, demonstrating that the proposed method allows the better reconstructing of initial energy density and target boundary.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use three image quality indexes: root mean square error (RMSE), structural similarity index (SSIM),^[66] and peak signal-to-noise ratio (PSNR) for quantitative analysis and the results are shown in Table 1. The proposed TAT-Net obtains lower RMSE and higher SSIM and PSNR than the DAS, which indicates that the TAT-Net has obvious advantages over the traditional method in reconstruction accuracy. Compared with the U-Net, the TAT-Net has a slight improvement in RMSE and PSNR, and a substantial improvement in SSIM, suggesting that the TAT-Net has better structural preservation. We also record the reconstruction times of different methods, and the average computing time of TAT-NET is about 8 times faster than that of DAS.

In the experiments, one to three saline-containing tubes are used as reconstruction target(s), and the TAT phantom dataset is obtained by using the TAT system shown in Fig. 2. The network is trained on the synthetic dataset and tested on the TAT phantom dataset. Although both datasets contain similar circular reconstruction targets, significant domain gap can be observed in the sinogram due to the TAT system factors. These experiments are used to confirm whether the TAT-Net has the ability to generalize the experimental dataset.

Figure 5 shows the reconstruction results of 1 to 3 circular targets from different methods. The reconstructed images by the DAS suffer severe artifacts and blur target edges. There also exists an annular reflection signal interference around the target. In the case of 3 targets, the larger target is slightly elliptical and splits into two parts in the DAS image, and the shapes of two lager targets also become elliptical in the TAT-Net image. This is, we believe, due to inhomogeneous microwave absorption caused by microwave line polarization, resulting in a split elliptical shape of the reconstructed target.^[67] Correspondingly, the sinogram of the larger target also splits into two lines. This type of data does not exist in our simulation training set, so the TAT-Net will experience distortion when reconstructing the data of such a type. A potential solution to this problem is to consider the signal change caused by microwave polarization in the training set simulation.

Zoom In Zoom Out Reset image size

The U-Net can reconstruct the position and shape of each target, but severe artifacts appear, especially in the 3-target case. By observation, it can be obviously seen that these artifacts and sinogram have very similar characteristics, and we suspect that the appearance of sinogram-like artifacts are due to the shortcut connection structure in the U-Net, and the shallow connections in U-Net can give rise to the decrease of the generalization ability and robustness of the model. We note that the TAT-Net removes the shortcut link structure compared with U-Net and provides the artifact-free images, which demonstrates our hypothesis.

The 1D line intensity profiles shown in Fig. 6 clearly demonstrate the stability of TAT-Net against disturbances from TAT system. In addition, the TAT-Net obtains a correct target shape based on real characteristics, which shows great robustness against domain gap.

Zoom In Zoom Out Reset image size

Since the training set only contains the data of 3 targets at most, we study the robustness of TAT-Net against more imaging targets. Figures 7 and 8 show the results on 4 and 5 circular targets in simulations and experiments, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In simulation, the TAT-Net can completely reconstruct 4 and 5 circular targets, however, some circles have slight elliptical distortion and shift in the 5 circular case. Unlike the distortion in Fig. 5, this problem caused mainly by the generalization error of the model. In the TAT phantom experiments, it can be seen that the DAS is susceptible to artifacts that produce a blurred image and U-Net is almost impossible to reconstruct correctly. By contrast, the TAT-Net can still give a satisfactory quality of image reconstruction, demonstrating its great potential in the more complex situations.

Most of the existing TAT reconstruction methods are qualitative, which means that the reconstructed image can only provide the structural characteristics and relative amplitude of the reconstruction target. The TAT-Net has the ability to be quantitatively reconstructed by taking the advantage of synthetic data, where the real initial energy density can be used as ground truth. In this subsection, we verify the ability to quantitatively reconstruct the energy loss density by using the TAT-Net. We study 64 sets of simulation samples and calculate the maximum value inside a reconstruction target as the initial energy loss density. Figure 9(a) shows the distributions of all the samples and the fitting line based on the least square method. It can be clearly seen that there is a very high correlation (R = 0.9858) between the TAT-Net and the ground truth; this demonstrates that the proposed method provides a stable and accurate quantitative energy density prediction.

Zoom In Zoom Out Reset image size

In the experiment, one saline-containing tube is placed at different positions and we can obtain a series of phantom data with different energy density distributions due to the inhomogeneous microwave energy field. The real energy density cannot be measured, and thus, we take the reconstructed images by DAS as reference and calculate the correlation (R = 0.9295) between TAT-Net and DAS. Figure 9(b) gives the results of 90 sets of experimental samples, which indirectly demonstrate the ability of the proposed method to implement the quantitative reconstruction.