A Comparative Study of Four Total Variational Regularization Reconstruction Algorithms for Sparse-View Photoacoustic Imaging

Photoacoustic imaging (PAI) is a new nonionizing, noninvasive biomedical imaging technology that has been employed to reconstruct the light absorption characteristics of biological tissues. The latest developments in compressed sensing (CS) technology have shown that it is possible to accurately reconstruct PAI images from sparse data, which can greatly reduce scanning time. This study focuses on the comparative analysis of different CS-based total variation regularization reconstruction algorithms, aimed at finding a method suitable for PAI image reconstruction. The performance of four total variation regularization algorithms is evaluated through the reconstruction experiment of sparse numerical simulation signal and agar phantom signal data. The evaluation parameters include the signal-to-noise ratio and normalized mean absolute error of the PAI image and the CPU time. The comparative results demonstrate that the TVAL3 algorithm can well balance the quality and efficiency of the reconstruction. The results of this study can provide some useful guidance for the development of the PAI sparse reconstruction algorithm.


Introduction
Photoacoustic imaging (PAI) is a novel noninvasive biomedical imaging modality with the capability of quantitatively imaging of light absorption characteristics of endogenous tissue chromophores that has grown tremendously in the last two decades [1][2][3][4]. As a hybrid imaging method, PAI combines strong optical contrast with high ultrasonic penetration [5][6][7]. And it has shown great potential in multiple clinical applications, including the breast imaging [8,9], dermatologic imaging [10,11], thyroid imaging [12,13], and imaging of the lymphatic system [14]. In PAI, the reconstruction of a high-quality image usually requires a large amount of signal data, which requires expensive electronic equipment or long data acquisition times. Moreover, in many clinical applications including ophthalmic imaging [15] and breast imaging [8], only incomplete data with limited angles can be accepted. Additionally, the conventional analytic methods usually reconstruct distorted images with strong artifacts when the data are insufficient or collected from few views. Therefore, the development and investigation of high-speed and high-quality PAI image reconstruction algorithms based on incomplete data is a popular research area of current interest [16][17][18]. To address the issue of insufficient information, the iterative algorithms for PAI reconstruction have been proposed to improve the quality of reconstructed images and reduce the time of data acquisition [19][20][21]. The incomplete data may arise from a variety of forms, but in this work, we focus on the sparse data problem in PAI with a circular measurement geometry.
Mathematically, the reconstruction of PAI images from sparse data can be regarded as a problem of solving underdetermined linear equations. By incorporating some prior information of the object or missing data, the iterative algorithms that can obtain more accurate reconstruction results at the cost of much more computing time have been developed for PAI [19][20][21]. One of them is based on the theory of compressed sensing (CS), which has attracted more and more attention due to its can recover sparse signals using much less measurements than advised by Shannon's sampling theory [22,23]. By using the L1magic convex optimization algorithm and sparse-view data, Provost and Lesage applied the CS theory to the PAI imaging for the first time [24]. The issue of loss of resolution and artifacts in the case of insufficient measurements can be solved by using random lighting via the SPGL1 algorithm [25]. Based on a set of highly sparse representations of noise signals, Bayesian CS theory was used to obtain PAI images [26]. The results of phantom and in vivo experiments showed that the CS method can effectively reduce the undersampling artifacts through the nonlinear conjugate gradient descent algorithm [27,28]. All these studies showed that the CS-based iterative reconstruction methods can significantly reduce the number of ultrasound transducers in the PAI imaging system and obtain high-quality reconstructed images with sparse data.
The total variation (TV) regularization and compressed sensing theory has a wide array of applications in biomedical imaging, for its good properties for preserving sharp edges and contours of objects. For example, a TV iterative shrinkage scheme is proposed to CS recovery of parallel MRI [29,30]. Jia et al. developed the GPU-based cone-beam CT reconstruction algorithm using noisy and reduced projection data via TV [31]. Liu et al. investigated the application of sparse Bayesian learning framework and in electrical impedance tomography [32,33]. This paper conducts a quantitative performance study on the application of CS-based TV regularization reconstruction algorithm in PAI. The total variation minimization by augmented Lagrangian and alternating direction algorithm (TVAL3) [34,35], the generic log-barrier algorithm (L1magic) [24,36], Nesterov's algorithm (NESTA) [37], and the two-step iterative shrinkage thresholding algorithm (TwIST) [38] are considered to solve the TV minimization problem. Based on sparsely sampled data, we evaluated the performance of the four reconstruction algorithms. Both the quality of the reconstructed images and the CPU runtime are investigated. The results of this study are expected to provide a suitable reconstruction algorithm for PAI that reduces the scanning time without reducing the quality of the reconstructed image.
The organization of this article is as follows. Section 2 briefly reviews the methods used in the study, including the photoacoustic theory, reconstruction algorithms, and evaluation criteria. The results of a comparative experimental study and discussion are presented in Section 3. In the last part of this article, we got some conclusions based on the results of numerical experiments and simulation experiments.

Photoacoustic Theory.
According to the theory of photoacoustic signal generation, the relationship between the ultrasonic pressure pðr, tÞ in a homogeneous medium and the absorption distribution AðrÞ obeys the following wave equation [1]: where r denotes the pixel coordinate, t means the time, c is the sound velocity, β is the thermal coefficient of volume expansion, C p is the specific heat, and δðtÞ is a delta function. By using the Green function, wave equation (1) can be solved, providing the forward problem [9].
where r 0 is the position of the ultrasound transducer. By taking the Fourier transform of equation (2) and denoting k = ω/c, where ω is the angular frequency, the forward problem in the temporal-frequency domain can be expressed as [24] During the experiment, we obtained the photoacoustic signal in the frequency domain by performing fast Fourier transform on the time domain signal. After discretisation, equation (3) can be represented as a linear equation: where Y ∈ R M denotes the vector for all pressure measurements in Fourier domain and X ∈ R N denotes the vector for the unknown reconstruction image. However, we can only observe inaccurate measurements Y = KX + ε, where ε ∈ R M denotes the modeling transducer noise. According to equation (3), the time-frequency domain measurement matrix can be written as where r m represents the transducer position, r ij indicates the coordinate of pixel, p is the transducer quantity, and q represents the number of sampling locations, respectively. The image reconstruction problem is to extract the absorption distribution X from the pressure data Y, which is usually ill-posed when there are fewer sampling points. Therefore, regularization techniques and prior information of the image are utilized to obtain a stable reconstruction process.

Total Variation
Regularization. So far, a variety of CS reconstruction algorithm has been developed, such as the greedy iterative algorithm, total variation (TV), and Bayesian framework. Among them, TV regularization has great success in image reconstruction due to its capability of keeping edges and boundaries. The discrete form of the isotropic TV for a grayscale image is the sum of the L 2 norm of the 2 Computational and Mathematical Methods in Medicine discrete gradient: where D i X = ðΔ h i X, Δ v i XÞ T is the discrete gradient of the gray image X and Δ h i and Δ v i denote the horizontal and vertical difference operators. The existing literature indicated that TV-based CS algorithms can accurately reconstruct images from sparsely sampled data [34][35][36][37][38]. The L 2 norm of D i X corresponds to the isotropic TV, or the L 1 norm of D i X corresponds to the anisotropic TV. In the field of PAI, the TV minimization optimization algorithm can reconstruct excellent image from few view data [39][40][41][42]. The noiseless discrete TV regularization model can be written as For the reconstruction with noise signals, we can solve the Rudin-Osher-Fatemi model as an alternative In this work, the TVAL3, L1magic, NESTA, and TwIST are considered to solve TV minimization problem (7) or (8). Comparative performance is assessed for both noiseless and noisy sparsely sampled data.

Evaluation Factors.
Three quantitative parameters are used to evaluate the performance of these four reconstruction algorithms from a sparse data set: the CPU time, the signal-to-noise ratio (SNR), and the normalized mean absolute error (NMAE). The CPU time was used to estimate the efficiency of the algorithm, and the SNR was applied to evaluate the quality of the reconstructed image, and the NMAE was used to quantify the reconstruction error between the gray image and its reconstruction. The SNR is defined by where X _ is the reconstruction and X is the mean intensity valve of X. The NMAE is defined as

Experiments and Results
In this article, we provide a simulation-based comparative performance study between these four TV regularization algorithms. The forward simulation and inverse reconstruction were conducted in 2D phantoms. The NCAT phantom and the blood vessel phantom are used in the comparisons to generate the photoacoustic signals. And the photoacoustic signals are generated by using equation (3). Figure 1(a) shows the NCAT phantom, and Figure 1(b) shows the blood vessel phantom. The size of the phantom is 42 mm × 42 mm with a resolution of 128 × 128 pixels. During the simulation, the sound speed is 1500 m/s and the diameter of the circular scan is 60 mm. To simulate the response of the ultrasonic transducer, at every sampling location, 64 randomly chosen k n /2πc's inside the ½0:1, 32 MHz window were employed to define the frequency domain projection matrix K. By rescaling the phantom gray values to ½0, 1, we acquired the frequency domain signals using the projection matrix.
The simulation experiments were carried out using MATLAB (MathWorks, Natick, MA) on a personal computer with an Intel core i7-4790 processor and 32 GB memory. The parameter ranges of the four algorithms are selected according to the literature 34 to 38, and the specific parameter values are manually adjusted. The same iteration stopping criteria kX k+1 − X k k/kX k k < 0:005 were used to be fair to compare the four algorithms. In the following section, the quality of the reconstructed images and quantitative comparisons are going to be discussed.
3.1. Sparse-View Reconstruction. The reconstructed images of the NCAT phantom using these four TV regularization algorithms and TV regularization model (7) from a set of the sparsely sampled signals, with 20, 30, and 40 positions, are shown in Figure 2. The first column displays the reconstructed images of the TVAL3 algorithm, and the second column show those of the L1magic algorithm. The third column is the reconstructed images of the NESTA algorithm, and the fourth column presents those of the TwIST algorithm. From Figure 2, it can be found that the image of the NCAT phantom has been reconstructed very well by the TVAL3 algorithm, even when the sampling number is reduced to 20. There are stripe and speckle noise in the image reconstructed by the L1magic algorithm, which affects the quality of the reconstruction. The image reconstructed by the NESTA algorithm is blurred, which affects the visual effect. Among the four algorithms, the TwIST algorithm performed the worst. Even if signals with 40-view angles are used, the TwIST algorithm still cannot obtain a better image. The TwIST algorithm still cannot obtain a better image under 40-view sampling circumstance. This experiment shows that the TVAL3 method is better than other three methods significantly in PAI image sparse reconstruction.
In order to better observe the differences between these four algorithms, we drew the pixel gray scales along the middle extraction lines of the reconstructed images. In Figure 3, the red dotted line and the black solid line are the pixel gray scales of the reconstructed image and the NCAT phantom, respectively. By observing the gray scales, we can get the similar conclusions mentioned above. Therefore, when these two evaluation factors are considered together, the TVAL3 algorithm achieves the best reconstruction performance, followed by the L1magic algorithm and the NESTA algorithm.
The quantitative evaluation parameters of the reconstructed images including the CPU time, the SNR, and the 3 Computational and Mathematical Methods in Medicine NMAE achieved by these four algorithms. From the CPU time in Table 1, we can know that there are large differences between the algorithm execution times: TwIST is about 1.3 times faster than NESTA, which itself can be roughly 2 times faster than TVAL3. And the L1magic algorithm has the longest running time. As the sampling number increases, the SNR data show an upward tendency and the NMAE values show a downward tendency. According to Table 1, we find that the TVAL3 algorithm provides the highest SNR value and the lowest NMAE value and the TwIST algorithm has the lowest SNR value and the highest NMAE value. NESTA and L1magic have the similar reconstruction performance. Considering the reconstruction efficiency and performance, the TVAL3 algorithm is the optimal algorithm for sparsely sampling PAI. Figure 4 illustrates the reconstructed images of the blood vessel phantom using these four TV regularization algorithms from a set of the sparsely sampled signals, with 20, 30, and 40 positions. Therein, the first column to the fourth column show the reconstruction results of the TVAL3 algorithm, the L1magic algorithm, the NESTA algorithm, and the TwIST algorithm, respectively. From Figure 4, we can see that the blood vessel phantom has been well reconstructed when using the TVAL3 algorithm and 40 position signals. The image quality of the reconstruction is reduced when 30 position signals are used. And the quality of image reconstructed by using 20 position signals becomes very poor. However, the reconstructed images contain noises and blurs when using the other three algorithms and 40 position signals, and the quality of the reconstructed images  This experiment fully proves that the TVAL3 method has the best reconstructed image quality. The detail distinguish ability can also be seen from the pixel gray scales along the middle extraction lines of the reconstructed images (the red line in Figure 1(b)). As can be seen from Figure 5, the TVAL3 algorithm has the best image detail reconstruction capability. The reconstructed image with L1magic algorithm contains more interference noise. The details of the image reconstructed by the NESTA algorithm are blurred. And the TwIST algorithm has the worst image detail reconstruction ability. According to We record the CPU time for these four algorithms in each experiment of the blood vessel phantom, as shown in Table 2. The CPU time of the TwIST algorithm is approximately 2.2 times longer than the NESTA algorithm, which itself is about 2.6 times faster than TVAL3. And the CPU time reconstructed by the L1magic algorithm is the longest. According to the SNR in Table 2, we can see that the TVAL3 algorithm can obtain a SNR of 30 dB using signals from 30views, while the SNRs of the reconstructed images of the other three algorithms are less than 30 dB using signals from  Considering the results of these two simulation experiments, we can draw a useful conclusion that the quality of reconstructed image with the TVAL3 algorithm is better than the other three algorithms. The artifacts and blurs emerge in the other three algorithms reconstructed images, and the quality of images is severely affected indicating that these three methods are not suitable for sparse reconstruction. From the simulation results, we can conclude that the TVAL3 algorithm is more suitable for reconstructing photoacoustic images under sparse sampling conditions than the other three algorithms.

3.2.
Robustness to the Noise. In the actual photoacoustic imaging process, the signal will be affected by Gaussian   Figure 5: Gray curves along the middle extraction line in Figure 1(b). The gray curves are reconstructed with these four algorithms from 20view, 30-view, and 40-view signals for the NCAT phantom. 6 Computational and Mathematical Methods in Medicine white noise from the system ultrasonic transducer and electronic devices. Therefore, it is very important for an algorithm to keep stable performance in the presence of noise. In order to analyze the robustness of those four methods, we added different levels of Gaussian noise to the signal. Results of noisy signal experiments are evaluated by SNR and NMAE. In Figure 6, the SNR and the NMAE indexes are represented as a function of the number of data. A good algorithm has a larger SNR and a smaller NMAE. Figure 6 shows the trend diagram of the SNR and NMAE between the NCAT phantom and the reconstruction obtained from those four methods and TV regularization model (8). We remind that higher SNR and lower NMAE indicate superior image quality. It is clear that the TVAL3 algorithm outperform the other three algorithms also in terms of the SNR and NMAE. Moreover, the TVAL3 algorithm can not only obtain better performance with fewer measurements but also improve SNRs faster than the other three algorithms.
It can be seen from Figure 7, when there is weak noise of 40 dB and 30 dB, those four methods achieve similar SNRs and NMAEs when using signals with less than 25 positions. Furthermore, the TVAL3 algorithm achieves the biggest SNR and the smallest NMAE with sampling locations of more than 40 and improves SNR faster than the other three algorithms. When there is strong noise of 20 dB and 10 dB,  7 Computational and Mathematical Methods in Medicine the NESTA algorithm and the TwIST algorithm have smaller SNRs and larger NMAEs, but the TVAL3 algorithm can still achieve good reconstruction results with sufficient measurements. In other words, the TVAL3 algorithm achieves good reconstruction of PAI images in noisy environments. Figure 8 shows the histogram of the influence of changing variance of Gaussian noise distribution on TVAL3. When the number of signals is small, there is not a great deal of difference in the SNR of reconstructed images with different levels of noise signals. While the number of signals exceeds 25, the SNR of the reconstructed images under different noise levels will vary greatly. It can also be seen from Figure 8 that the SNR of the reconstructed image can be improved by increasing the number of signals when the noise level is low, and the SNR of the reconstructed image is low when the noise level is high. Therefore, the biggest challenge in practical experiments is how to minimize the impact of noise. Signal averaging technology is needed to eliminate noise interference and obtain more reliable and effective signal data.

In Vitro Experiments.
We also compared those four reconstruction algorithms through in vitro experiment. Figure 9(a) shows the schematic diagram of the circular scanning experimental system, which is modified from our previous article. A Q-switched 532 nm Nd : YAG laser was used as the source of light. The single-element ultrasound transducer (V309 Panametrics) with a central frequency of 5 MHz and a diameter of 12.7 mm was used to receive the ultrasound signals. At each signal acquisition position, the photoacoustic signals were first amplified by a pulse  Computational and Mathematical Methods in Medicine amplifier, then recorded with an oscilloscope (MSO4000B, Tektronix), and finally input into a personal computer for signal processing and image reconstruction. The imaged phantom we used in the experiment is made by gelatin cylinder. The agar phantom was made by mixing 1% lipid, 6% gelatin, and 93% water to simulate biological tissue. Figure 9(b) is the photograph of the phantom. The diameter of the phantom is 20 mm. One graphite rod with a diameter of 0.5 mm and two hairs of lengths of 4 mm were embedded as the optical absorbers. In the phantom experiment, 60-view data are selected for reconstruction.
In the experiment, 40-view and 80-view signals that are uniformly distributed at a 360°curve are collected to reconstruct images. Figure 10 displays the images that were reconstructed from the sampling data sets utilizing TVAL3, L1magic, NESTA, and TwIST, respectively. The second row of Figure 10 is reconstructed from 80-view data. When the sampling signals is sufficient, the TVAL3, L1magic, and NESTA methods can construct good-quality images, and the TVAL3 method has the best reconstructed image quality. However, the quality of the reconstructed image with the TwIST method is still relatively poor. When we reconstruct the image with 40-view data (first row of Figure 8), the reconstructions are seriously affected and there are more noise in the images. Only the image reconstructed by the TVAL3 algorithm is relatively clear, while the noises in the reconstructions of the L1magic and NESTA methods affect the identification of the boundary of the phantom.

Conclusion
In this paper, we have tested and evaluated four reconstruction algorithms (TVAL3, L1magic, NESTA and TwIST) for PAI images with sparsely sampled data. The numerical simulations demonstrate that the TVAL3 method has the most accurate reconstruction performance under sparse sampling while the TwIST method has the worst detail reconstruction ability. In terms of CPU time, the TwIST method needs the least amount of CPU time, and the L1magic method requires the longest CPU time. In terms of the reconstructed image quality, the TVAL3 algorithm has better accuracy and antinoise ability than the other three algorithms. In conclusion, the TVAL3 algorithm has the best image quality and requires less CPU time, which provides a good balance between the accuracy and efficiency of the reconstruction. We believe that the findings of this research will provide

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.