Abstract
Photoacoustic imaging involves reconstructing an estimation of the absorbed energy density distribution from measured ultrasound data. The reconstruction task based on incomplete and noisy experimental data is usually an ill-posed problem that requires regularization to obtain meaningful solutions. The purpose of the work is to propose an elastic network (EN) model to improve the quality of reconstructed photoacoustic images. To evaluate the performance of the proposed method, a series of numerical simulations and tissue-mimicking phantom experiments are performed. The experiment results indicate that, compared with the -norm and -normbased regularization methods with different numerical phantoms, Gaussian noise of 10-50 dB, and different regularization parameters, the EN method with has better image quality, calculation speed, and antinoise ability.
1. Introduction
During the recent 20 years, photoacoustic imaging (PAI) has received wide spread attention because it combines the advantages of both optical and ultrasound imaging [1–3]. PAI can obtain the quantitative information of the light absorption coefficient of biological tissue and hence has been successfully used in various clinical applications, including breast imaging [4, 5], skin imaging [6, 7], cardiovascular imaging [8], and endoscopy [9, 10]. Photoacoustic effect is excited when a short pulse laser irradiates the biological tissue. And the optically absorbing targets in the tissue will produce the photoacoustic signals, which can be received by the ultrasound transducers placed around the tissue. These acquired signals can be utilized to calculate the spatial light absorption distribution inside the tissue through a certain image reconstruction method [11, 12].
The reconstruction algorithm has an important effect on image quality. Traditional reconstruction algorithms, such as FBP and delay and sum algorithms [11–13], are widely used because of their accuracy and convenience. However, these algorithms need to collect complete data to reconstruct a high-quality image. Furthermore, only a limited range of photoacoustic signals can only be obtained in many experiments [14–16]. In this case, the PAI images reconstructed by the traditional algorithms usually suffer from streaking artifacts and edge blurring. Generally, PAI reconstruction with incomplete data is an ill-posed problem, and regularization techniques must be applied to suppress noise and artifacts in the reconstructed images. By adding some prior knowledge or reasonable regularization constraints, the model-based iterative reconstruction algorithms that can further improve the image quality have been developed for PAI [17–19]. One of the algorithms is based on -norm regularization, where the -norm of the light absorption distribution is applied to constrain the PAI reconstruction. The most popular -norm regularization approach that can achieve better prediction performance has been used in PAI [20, 21]. However, -norm regularization often generates oversmooth solutions. As an alternative, -norm regularization attracted enormous research interest in PAI due to its sparsity-inducing property [22–24]. Nevertheless, the -norm regularization cannot reconstruct an image with the least signals when applied to compress sensing. To be able to achieve more sparse results, the nonconvex -norm regularizations for PAI have been investigated, which provides better reconstruction results over -norm regularizations [25–27]. The above studies have confirmed that -norm regularization has been widely applied in PAI. However, the -norm regularization is a nonconvex optimization problem, which is difficult to solve quickly and efficiently.
Recently, the elastic net (EN) regularization is developed for solving the ill-posed inverse problem of -norm regularization [28]. As a combination and compromise between - and -norm regularizations, the EN method is useful when the total number of image pixels is much larger than the number of observed signals, which has been applied to many biomedical imaging fields. For example, Majumdar and Ward have used EN to reconstruct MRI images based on the undersampling sparse -space data [29]. Wang et al. applied the adaptive parameter search EN algorithm to calculate the inverse problem of the fluorescence molecular tomography [30]. Causin et al. investigated the application of EN regularization technique in diffuse optical tomography [31]. All the above papers have verified that EN regularization technique can produce a sparse model with good prediction accuracy.
In this paper, the EN regularization is used to optimize the PAI image reconstruction problem. The paper is organized as follows. The photoacoustic theory and the reconstruction algorithm based on the EN regularization are presented in the second section. The third section introduces the experimental processes and results of the numerical simulation and the tissue-mimicking phantom. In Section 4, some conclusions are given.
2. Theory and Method
2.1. Photoacoustic Theory
According to the photoacoustic signal generation theory, the relationship between pressure at position and time and the laser energy deposition in a homogeneous acoustic medium can be described by the following: where is the acoustic speed in the medium; and denote the specific heat and the thermal coefficient of volume expansion; and represents the shape of the laser pulse that can be expressed as a delta function. Based on Green’s function, the solution to the time domain forward problem can be expressed as [1] where is the location of the ultrasonic detector. By denoting and taking the Fourier transform of time , it is easy to obtain the frequency spectrum of the photoacoustic pressure
The purpose of the inverse PAI problem is to obtain an approximation of the absorption distribution from a set of sampled data . The forward problem can be described using the following matrix relation: where is a column vector representing the acoustic pressure and is the number of measurements; denotes the vector of the unknown reconstruction image and denotes the number of pixels; and represents the noise. is the projection matrix in the temporal frequency domain with the size of that can be can be discretized as where is the position of the ultrasonic transducer and is the total number of location, represents the frequency and is the number of samples, and indicates the image pixel coordinates, respectively. Generally speaking, the well-known FBP algorithm can obtain better reconstruction accuracy when the signal is sufficient. However, the FBP algorithm has poor prediction accuracy when the data is incomplete, and the regularization methods are needed to improve the reconstruction quality.
2.2. -Norm Minimization-Based Iterative Reconstruction Algorithm
The PAI image reconstruction with insufficient measurements is essentially an ill-posed inverse source problem. It is difficult for the analytic reconstruction algorithm to reconstruct high-quality images in limited data settings. Therefore, regularization is employed to suppress noise and artifacts and yields an acceptable image. The -norm-based Tikhonov regularization method is often employed to solve the ill-posed problem [21, 32]. We can get the regularized solution of Equation (4) in the framework of Tikhonov by minimizing the following functions: where is a regularization parameter. And the regularization operator is usually chosen to be the identity matrix or finite difference operators. The Tikhonov regularization can be solved effectively by LSQR algorithm. However, the oversmoothness of the Tikhonov regularization solution will lead to the loss of details in the reconstructed image.
The compressed sensing (CS) theory suggests that an image can be reconstructed exactly from insufficient measurements if it is sparse or can be sparsely represented in an appropriate basis. Fortunately, most of the medical image can be sparsely represented under a suitable sparse transform basis , in which is the sparse transform coefficient. It has been proved that the PAI image is sparse in a discrete wavelet basis and numerical derivative basis [22]. In this paper, is defined as four-level symmetric wavelet transform. In order to reconstruct PAI images with insufficient measurements by using the theory of compressed sensing, the following -norm optimization problem can be solved:
At present, the most common solution of Equation (7) is to employ the CS techniques, such as L1-MAGIC [22], YALL1 [23], and SPGL1 [24]. However, -norm-based LASSO regularization usually has some serious flaws such as oversparseness. In recent years, Zou and Hastie have proposed the elastic network regularization methods to improve the problem that the traditional regularization methods are too smooth or too sparse [28].
The functional interpretation of the solution with elastic net regularization term is shown below. Aiming at the problem that the traditional regularization method is oversmooth or oversparse, the elastic network regularization method was proposed. The framework is as follows: where is the overall regularization parameter, denotes the convex combination weight, which decides the weights of the - and -norm terms. is the EN penalty, which is a compromise between LASSO and the Tikhonov regularization. This penalty possesses the advantages of both regularizations and meets properly the requirement of image sparsity and smoothness. If , Equation (8) becomes the -norm-based LASSO regularization; If , Equation (8) will become the -norm-based Tikhonov regularization. The quality of image reconstruction by the EN method is usually related to the selection of regularization parameters. It is very important to choose the values of parameters and correctly.
2.3. The Evalution Factors
To quantitatively evaluate the efficiency and accuracy of PAI reconstruction, the CPU running time, the normalized mean absolute error (NMAE), and the peak signal-to-noise ratio (PSNR) were used as quantitative factors. The CPU runtime is used to evaluate the computational efficiency of the PAI reconstruction algorithm. The PSNR is applied to estimate the image quality, and the NMAE is employed to quantify the reconstruction error. The PSNR can be defined by where is the reconstruction image and means the maximum value of pixels, which in our simulation is 1. The NMAE is defined as
3. Experiment and Result
In this section, multiple numerical simulations and applications were performed to validate the effectiveness of the EN method. And the -norm-based Tikhonov regularization method [33] and the -norm-based regularization method SPGL1 [34] are used to compare with the EN method, where the Tikhonov regularization method is a special case of the EN method corresponding to . Both forward projection and backward reconstruction are performed in 2D, where the object to be imaged is approximately contained in a thin plate. The sparse transform operator is set to four-level symmetric wavelet transform using the Rice Wavelet Toolbox. All the MATLAB programs are performed on a desktop computer with a 3.6 GHz CPU and 32 GB memory. In this paper, we use the MATLAB package glmnet to select the appropriate parameters of the EN [35].
3.1. Reconstruction from Simulated Sparse-View Data
Figure 1 shows the breast phantom and the blood vessel phantom. The phantoms and Equation (4) are used to produce the photoacoustic signals. The size of the phantom is and the resolution is pixels. During the simulation experiment, the diameter of the circular scan of the ultrasonic transducer is 45 mm, and the sound velocity of the ultrasonic is 1500 m/s. At each sampling position of the ultrasound detector, 64 samples randomly selected in the (0.15, 4) MHz window are used to define the projection matrix using Equation (4). By normalizing the gray value of the phantom to (0, 1), the simulation signal is obtained by using the frequency domain projection matrix. The frequency domain measurement data can be generated by using .
Figure 2 shows the experimental results of breast phantom using these three algorithms. It can be seen that all three methods achieve good reconstruction results when using 60 position signals. Moreover, the reconstruction capability of the EN method is better than that of the Tikhonov and SPGL1 methods in the visual sense. When using 40 sampled signals, the reconstructed image of these three methods contain large noise and artifacts. In addition, Tikhonov has the worst image quality indices, followed by SPGL1. The quality of photoacoustic images reconstructed by all these three methods is poor when 20 sampled signals are used. As can be seen from Figure 2, the EN method can obtain more accurate images when the sampling number is sufficient.
In order to verify the general applicability of the EN method, we choose the blood vessel phantom as the initial optical deposition to additionally compare these three algorithms. And this experiment has the same simulation environment as the breast phantom experiment. As can be seen from Figure 3, Tikhonov has darker and noisy reconstructed images when small amounts of signals are used. In addition, the SPGL1 method is unable to reconstruct high-quality photoacoustic images based on signals from less than 30 locations, while the EN method can reconstruct clear photoacoustic images using signals from 20 sampling locations. From the above two simulation results, it can be seen that the EN method can effectively remove noises and preserve edges.
As can be seen from the first column of Figure 4, the CPU time becomes larger when the number of measurements becomes larger. The running time of two phantom experiments is similar. The Tikhonov method takes the longest time to reconstruct, and the EN method has the shortest reconstruction time. From the second column of Figure 4, we can see that the larger the number of samples, the lower the NMAE. In the breast phantom experiment, 50-60 sampled signals are required to obtain a small NMAE value, while in the blood vessel phantom experiment, only 30 sampled signals are required to obtain a satisfactory NMAE value. It can be seen from Figure 4 that the Tikhonov method fail to generate smaller NMAE and larger PSNR even with 60 sampled signals. In the blood vessel phantom experiment, when the number of sampled signals is sufficient, the NMAE and PSNR obtained by the SPGL1 method are similar to those obtained by the EN method. However, when fewer signals are used, the PSNR of the SPGL1 method is much lower than that of the EN method. The above two simulations verify that the EN with and yields a higher time resolution and a higher reconstruction quality than the Tikhonov and SPGL1 methods.
3.2. Antinoise Ability Experiment
Noise is likely to be added during photoacoustic signal acquisition. To assess the accuracy and stability of the EN methods, we added various levels of noise to the simulated data and investigate the effects of noise on the reconstructed images. A stable algorithm has higher PSNRs and smaller NMAEs. The numerical results of the PSNR and NMAE from the blood vessel phantom experiment are shown in Figure 5. As can be seen from Figure 5, the EN method achieves maximum PSNR and minimum NMAE at various noise levels. The PSNRs of the Tikhonov methods lightly drop with the increase of noise level, and the PSNRs of theSPGL1 algorithm decrease the most, while the PSNRs of the EN algorithm are essentially unchanged. The results of the antinoise ability experiments confirmed that EN is the most robust of these approaches. And EN with has better noise robustness than EN with .
3.3. Parameter Investigation
In this subsection, we investigate the impact of the regularization parameter on image quality reconstructed by the EN method. Here we select 30-view simulated data from the blood vessel phantom to further investigate the parameter settings. And the value of ranges from 0 to 1. We calculate the PSNR value of the reconstructed image, and the line charts of PSNRs are shown in Figure 6. The numerical simulation experiment shows that, for noiseless data case, the PSNR value becomes larger and larger when the alpha value changes from 0 to 1, and it changes very little when the alpha value is from 0.4 to 1. Nevertheless, for noisy data case, EN with has the maximum PSNR value, that is, has better noise robustness.
(a)
(b)
3.4. Tissue-Mimicking Phantom Experiment
We performed the tissue-mimicking phantom experiment to evaluate the practicability of the EN approach. In Figure 7(a), the schematic of experimental set-up is shown. A Q-switched 532 nm Nd:YAG laser with a frequency resolution of 10 Hz was applied as the light source. The input laser pulse was amplified by a concave lens, homogenized with ground glass, and then irradiated onto a sample made of agar and black carbon sticks. A 5 MHz single-element ultrasonic transducer with a diameter of 12.7 mm (V309, Panametrics) was used to receive the photoacoustic signal. The transducer and the sample are submerged in a water tank to couple the photoacoustic waves to the transducer. A stepper motor (PMC100-3) controls the transducer to rotate around the sample for sampling. The rotation radius of the transducer is 40 mm. At every sampling point, the ultrasound signal was first amplified by a Panametrics pulse amplifier and then captured and averaged 30 times by an oscilloscope (MSO4000B; Tektronix). A personal computer is used to control the stepper motors and signal acquisition [36].
(a)
(b)
Figure 7(b) is the cross-sectional view of a cylinder agar phantom containing two carbon rod absorbers. The radius of the sample is 10 mm. Two carbon rods with a 0.5 mm diameter and lengths of 5 mm and 10 mm were embedded in the phantom as the optical absorbers. In the reconstruction experiment, signals at 40 and 80 positions evenly distributed on the circumference are used. Figure 8 was reconstructed by 64 frequency samples randomly chosen inside the (0.25, 5) MHz window. And the transducer response was restricted to certain value 1.
Figure 8 shows the results of agar phantom experiment. The images are constructed by the Tikhonov, SPGL1, and EN methods, respectively. There are too many noises in the images reconstructed by the Tikhonov and SPGL1 methods using 40-view signals, which cause all the details of the phantom to be suppressed. Although the image reconstructed by the EN method with 40-view data also has a lot of noise, the approximate outline of the phantom can be seen. The second row of Figure 8 is reconstructed from 80-view signals. When the number of sample points is sufficient, all three methods are feasible. The EN method with achieves the best results in terms of reconstruction accuracy, visual effects, and noise robustness among the three algorithms.
4. Conclusion
The principal purpose of this paper was to evaluate the application of the EN method in PAI. Based on CS theory, the EN method can reconstruct photoacoustic images using a small amount of data. To evaluate the reconstruction performance of EN method, we compare it with the Tikhonov and the SPGL1 in terms of visualization and performance indicators. The reconstruction results show that the EN method provides good imaging quality; it also simultaneously has acceptable time efficiency and better robustness. And the EN with has better noise robustness than EN with . Future work will focus on the biomedical applications of the EN method.
Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the Natural Science Foundation of Shandong Province (ZR2020QF025, ZR202102270451), the National Natural Science Foundation of China (61976110, 62176112), and the Discipline with Strong Characteristics of Liaocheng University—Intelligent Science and Technology under Grants 319462208.