A Photoacoustic Imaging Algorithm Based on Regularized Smoothed L0 Norm Minimization

The recently emerging technique of sparse reconstruction has received much attention in the field of photoacoustic imaging (PAI). Compressed sensing (CS) has large potential in efficiently reconstructing high-quality PAI images with sparse sampling signal. In this article, we propose a CS-based error-tolerant regularized smooth L0 (ReSL0) algorithm for PAI image reconstruction, which has the same computational advantages as the SL0 algorithm while having a higher degree of immunity to inaccuracy caused by noise. In order to evaluate the performance of the ReSL0 algorithm, we reconstruct the simulated dataset obtained from three phantoms. In addition, a real experimental dataset from agar phantom is also used to verify the effectiveness of the ReSL0 algorithm. Compared to three L0 norm, L1 norm, and TV norm-based CS algorithms for signal recovery and image reconstruction, experiments demonstrated that the ReSL0 algorithm provides a good balance between the quality and efficiency of reconstructions. Furthermore, the PSNR of the reconstructed image calculated by the introduced method was better than the other three methods. In particular, it can notably improve reconstruction quality in the case of noisy measurement.


Introduction
Photoacoustic imaging (PAI) is a new noninvasive and nonionizing biomedical imaging method which gains rapid development in the last two decades [1][2][3][4]. As a hybrid technique, PAI attains the advantages of the contrast of optical imaging and high resolution of ultrasonic imaging. In particular, PAI has shown great potential in many applications, including living mouse brain vasculature imaging [3,5,6], human skin imaging [7,8], and the treatment and diagnosis of cancer [9][10][11]. When the biological tissue is irradiated by nanosecond laser pulses, the photoacoustic signals will be generated in the tissue's light absorption region and measured by ultrasonic transducers. Afterwards, the distribution of laser energy absorption can be calculated from the photoacoustic signals by using the analytical or iterative algorithms.
The reconstruction algorithm is an important factor affecting the quality of PAI imaging, and the accurate and efficient reconstruction algorithms are of great significance. Analytic algorithms such as the filtered back-projection algorithm [12][13][14], the deconvolution reconstruction algorithm [15], and the time-reversal imaging algorithm [16] have already been used extensively owing to their accuracy and implementation convenience. However, the analytic algorithms can only reconstruct an accurate image with complete measurement data from all directions, which requires a long scan time or complex electronic equipment. When the data are insufficient, the analytic algorithms are no longer effective. In many PAI applications, such as breast imaging and ophthalmic imaging, only incomplete data can be accepted. There has not been accurate PAI reconstruction formula with incomplete data yet. Therefore, the development of high-speed and highquality PAI image reconstruction algorithms based on incomplete data is an active research topic recently. The incomplete data may arise from various forms, but in this work, we only consider the sparse-view PAI reconstruction problem.
Mathematically, image reconstruction with sparse-view incomplete data can be thought of as an underdetermined linear system. By making some constraints, the iterative reconstruction algorithm that gives a more accurate result at the expense of much more computing time has been developed for sparse-view PAI [17][18][19]. One type of them is based on the compressive sensing (CS) theory, which has drawn increasing attention due to its can recover sparse signals under sampling rate far lower than the Nyquist rate [20,21]. By using the L1magic convex optimization algorithm and sparse-view data, Provost and Lesage applied CS to PAI [22,23]. The problem of artifacts and loss of resolution in the case of insufficient measurements can be solved by using patterned excitation light via the SPGL1 algorithm [24]. Bayesian CS method was used to simplify the PAI system [25]. The experimental results indicated that the CS-based reconstruction method can reduce undersampling artifacts effectively through the nonlinear conjugate gradient descent algorithm [26]. The above studies indicated that the CS methods can reduce the number of ultrasonic transducers and get high-quality reconstruction results with sparse-view data.
To date, most CS applications within sparse-view PAI imaging have centered on the L 1 norm minimization problem [22,24,27] and the total variation (TV) minimization problem [28][29][30][31]. Encouraged by the advantages in edge preservation and finer structure recovering, the L 0 regularizations were introduced to computed tomography (CT) [32][33][34][35], magnetic resonance imaging (MRI) [36], and PAI [37]. Bioucas-Dias and Figueiredo proposed a smoothed L 0 norm (SL0) [38] algorithm that directly minimizes the L 0 norm, which combines the advantages of the high precision of convex optimization and the rapidity of greedy algorithm. Mozaffarzadeh et al. have shown that the SL0 algorithm can provide a higher PAI image quality while a low number of transducers were [37]. However, the SL0 algorithm does not consider the noise in the measurement data. Based on the SL0 algorithm, Bu et al. have adopted a regularization term to tolerate errors which can achieve the same computational efficiency as the SL0 algorithm while having better noise robustness. This inspired us to apply the regularized smooth L 0 (ReSL0) algorithm to reconstruct the sparseview PAI image. In this paper, we studied the use of the ReSL0 for L 0 norm minimization CS problems caused by sparse-view PAI reconstruction. The proposed algorithm is able to provide more accurate result under the sparse-view situation. For the purpose of verifying the capability of the ReSL0 algorithm, the ReSL0 algorithm compares with three L 0 norm, L 1 norm, and TV norm-based CS algorithms for signal recovery and image reconstruction.

Photoacoustic Theory.
Based on the photoacoustic signal generation theory, the relationship between the acoustic pressure pðr, tÞ and the absorbed energy density AðrÞ obeys the following wave equation [39].
With Fourier transform, the forward projection problem can be expressed in the time-frequency domain as During the experiment, the frequency domain data can be obtained by applying fast Fourier transform to the time domain measurements. To numerically model the above forward problem, the unknown reconstruction image AðrÞ is reshaped into one-dimensional long vector X ∈ R N , and a column vector Y ∈ R M can be used to represent all the temporal-frequency pressure pðr 0 , kÞ. According to Eq. (3), the temporal-frequency domain measurement matrix can be designed as [22].
which is determined by the geometry and the grid shape of the unknown image. In Eq. (4), r m indicates the position of the transducer, r ij denotes the pixel coordinates of the image, p is the number of ultrasonic sensors, and q indicates the number of sampling points, respectively. Then, Eq. (3) can be expressed as Y = KX. In PAI, the goal of image reconstruction is to reconstruct X through the pressure data Y.

Compress Sensing Application in PAI.
According to the theory of compress sensing, an image can be reconstructed when it or its transformation is sparse. Fortunately, most medical images can be considered sparse with a sparse transform basis Ψ : X = Ψθ, where θ ∈ R N indicates the sparse coefficient. Provost and Lesage have shown that there is a sparse transform basis in which the PAI image is compressible [22]. By denoting A = KΨ, the problem of PAI image reconstruction can be solved by the following optimization model.
where EðXÞ is the regularization function. Accroding to the literature [22], the matrix A obtained from the product between the forward operator K in the Fourier domain with a wavelets basis Ψ will be a CS-matrix. And the wavelet basis showed the best properties among generic bases that can represent sparsely images obtained via PAI. We note that most of the current CS-based reconstruction algorithms explore the prior knowledge that the PAI image is sparse or sparse in the transform domains. And the regularization term in Eq. (5) usually denoted by the L 1 norm of the image's sparsity transform coefficient, the total variation norm of the image, and so forth. For example, the constrained L 1 norm minimization can be applied to reconstruct the PAI image.

Molecular Imaging
PAI, the TV minimization optimization algorithm can reconstruct excellent image from few-view data [28][29][30][31]. In this work, the two-step iterative shrinkage thresholding (TwIST) algorithm is considered for image reconstruction [31,38]. where kθk 0 is the L 0 norm of θ [37,41]. In this algorithm, a continuous function was used to approximate kθk 0 instead of minimize the L 0 norm directly, which can be written as follows: f σ ðθÞ = exp ð−θ 2 /σ 2 Þ. It should have a parameter σ which determines the quality of the approximation. Considering Recently, Mozaffarzadeh et al. illustrate that SL0 offers higher-quality PAI images in comparison with L 1 normbased basis pursuit method while a low number of transducers were [37]. However, we can only observe inaccurate measurements Y = KX + e, which means there are errors between Y and KX. And the property of the SL0 decreases significantly due to the equality constraint Y = KX. In order to resolve this problem, we added observation noise to the forward model.
where e ∈ R M denotes a vector indicating the modeling transducer noise. And a regularized SL0 (ReSL0) method was used to solve the above problem mentioned in Eq. (10) [42,43], which can be written as follows: ReSL0 transforms the equality constraint in Eq. (8) into inequality constraint allowing certain error tolerance. The ReSL0 algorithm includes two nested iterations. The initial value σ 0 is set to 4 times of the maximum absolute value of the sparse coefficients in the outer loop, and the next values σ j = ρσ j−1 , j ≥ 1, where ρ is chosen between 0.5 and 1 in the experiment [41]. For every σ, the internal circle is responsible for finding the maximum value of F σ ðθÞ on set fθjkY-Aθk 2 < δg. The internal loop consists of iterations of the form θ ← θ + μσ 2 ∇F σ ðθÞ, followed by solving the optimization problem: Using Lagrange multipliers, this minimization results in [42].

Experiment and Result
Although the existing CS-based PAI reconstruction algorithms provide better results, the accuracy and efficiency of sparse-view PAI reconstruction still need further improvement. In this section, we provided a variety of simulations and in vitro applications to illustrate the advantages and efficiency of the ReSL0 method for sparse-view PAI reconstruction. The forward simulation and inverse reconstruction were conducted in 2D phantoms and images. 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. The simulation experiments were carried out using Matlab (version 7.8) on a PC with a 8 GHz CPU and 32 GB memory.
By using the wavelet transform and different sparse regularization methods, the numerical simulations have been carried out on the Sheep-Logan phantom, the blood vessel phantom, and the standard General Electric resolution phantom ( Figure 1) with a 128 × 128 resolution corresponding simulation area of 30 mm × 30 mm. During the simulation, the sound speed was 1500 m/s. A single ultrasonic transducer was used to receive the signals. To simulate the frequency response of the transducer, at every detection position, 128 randomly selected in the window [0.2, 3] MHz were used to define the projection matrix K. By adjusting the phantom gray values to ½0, 1, we obtained the ultrasonic pressure Y using the projection matrix K.

Reconstruction from Simulated Sparse-View Data.
We compared ReSL0 for the solution model (10) with SPGL1 for the solution model (6), TwIST for the solution model (7), and SL0 for the solution model (8). Figure 2 demonstrates the reconstruction results of the Sheep-Logan phantom by using these four methods. From the first row of Figure 2, we can see that the images of these four methods are strongly affected when 20 position samples are used. When the number of sampling is 35, the images reconstructed by the TwIST and SL0 methods contain artifacts and distortions, and the quality of the reconstructed images by the RESL0 and SPGL1 methods is better than these two methods. As the signal acquisition position reaches 50, the image reconstructed by the TwIST method still contains many noises, and the other three methods can reconstruct high-quality images. So we can conclude that the L 0 norm and L 1 norm-based CS algorithms can obtain more accurate image with sparse-view signals.
The quantitative evaluation of the reconstructed images including the CPU times, the peak signal-to-noise ratio (PSNR), and the normalized mean absolute error (NMAE) achieved by each of the algorithms is presented in Table 1. The PSNR (in dB) is defined as 10 log 10 ðN ⋅ MAX 2 I / kX −Xk 2 2 Þ and the NMAE is defined as 100 × kX − Xk 2 / kXk 2 , where MAX I means the maximum pixel value of the image andX denotes the estimate of the original. The CPU running time was used to estimate time complexity. The PSNR and the NMAE were used for quality evaluation of the reconstructed image.
According to Table 1, we can infer that there are great differences between the CPU times: TwIST is about 1.5 times faster than ReSL0, which itself can be roughly 7 times faster than SPGL1. And the SL0 has the same computational complexity as the ReSL0. From Table 1, we can see that those four methods achieve similar PSNRs and NMAEs when using signals at less than 35 locations. Table 1 also shows that the ReSL0 outperforms other three methods when using signals with more than 35 locations. Moreover, if we consider a method cannot reconstruct a high-quality image as the PSNR is less than 30, then we can conclude that the TwIST fails for all experiments, the SL0 fails when the number of sampling positions is less than 45, while the SPGL1 and ReSL0 algorithms fail when the number of sampling positions is less than 40. The numerical results in Table 1 indicate that the ReSL0 algorithm can reconstruct high-quality images with Considering the universality of the ReSL0 method, the blood vessel phantom was used as the initial energy density to additionally compare these four algorithms. The conditions of this experiment were the same as the Sheep-Logan experiment. Figure 3 shows the reconstruction images of the blood vessel phantom with these four algorithms. It can be seen from Figure 3 that all these four algorithms can reconstruct accurate images by using photoacoustic signals from 40 sampling positions. When the signal of 30 sampling positions is used, all reconstructed images contain much noise, which has a certain influence on the identification of the blood vessel structure. Among these four algorithms, the ReSL0 algorithm has the best reconstruction quality, which can be observed by the middle extraction lines from the reconstructed images. However, when the sampling position is 20, none of the four algorithms can get good reconstructed image. Table 2 shows the numerical results of the blood vessel for those four methods. Since the reconstructed image sizes of the two experiments are the same, the running time of    Table 2, we can see that the TwIST method can obtain a larger PSNR and a smaller NMAE when the number of samples is less than 30. However, the NMAE of the TwIST method is only slightly improved when the sampling position is greater than 30. And the ReSL0 method can achieve the largest PSNR and the smallest NMAE when the number of samples is greater than 30. According to Table 2, we can conclude that the TwIST method has better accuracy and efficiency when the signal number is small and the ReSL0 method obtains better quality PAI images with sufficient measurements.
In order to further verify the effectiveness of the RESL0 algorithm, a more complex and challenging standard General Electric resolution phantom is also used for simulation. The reconstruction results of the four algorithms are shown in Figure 4. In the first column of Figure 4, the images reconstructed by TwIST contain severe noises, which seriously affect the image quality. The second and third columns of Figure 4 are reconstructed by SPGL1 and SL0, respectively. It can be observed that both methods can suppress noise, but the quality of the reconstructed image can be further improved. The last column of Figure 4 shows the results of ReSL0, which obtain the best reconstruction quality among the four methods. Table 3 shows the numerical results of the standard General Electric resolution phantom for those four methods. As can be seen from Table 3, the TwIST algorithm has the shortest CPU time and the CPU time of the ReSL0 algorithm is slightly slower than the SL0 algorithm, while the CPU time of the SPGL1 algorithm is always much bigger than the other three algorithms. The TwIST algorithm has the worst PSNR and the biggest NMAE, and the PSNR value of images reconstructed by the ReSL0 algorithm is the highest, and the NMAE value of images reconstructed by the ReSL0 algorithm is the lowest. From the above three experiments, we can arrive at the conclusion that the ReSL0 algorithm outperforms the other three algorithms under sparse sampling condition.

Comparison for Different
Variances of the Noise. In order to quantify the influence of noise, the PSNR and NMAE of  10 Molecular Imaging the reconstruction images of the Sheep-Logan phantom experiment are calculated and shown in Figure 5. A good algorithm has a larger PSNR and a smaller NMAE. As can be seen from Figure 5, when there is weak noise of 50 dB and 40 dB, those four methods achieve similar PSNRs and NMAEs when using signals with less than 35 positions. Furthermore, the ReSL0 algorithm achieves the biggest PSNR and the smallest NMAE with sampling locations of more than 40 and improves PSNR faster than the other three algorithms. When there is strong noise of 30 dB and 20 dB, the SL0 algorithm and the TwIST algorithm have smaller PSNRs and larger NMAEs, but the ReSL0 algorithm can still achieve good reconstruction results with sufficient measurements. In other words, the ReSL0 algorithm achieves good reconstruction of PAI images in noisy environments. Figure 6 provides the PSNR and NMAE tendency chart of the reconstruction results of the blood vessel phantom experiment with noisy observation. When the number of sampling locations is greater than 25, the ReSL0 method can achieve the largest PSNR and the smallest NMAE. These results are consistent with the results of the noiseless environment as shown in Table 2. However, the performance of the SL0 algorithm degrades rapidly due to the equality constraint in Eq. (8). According to the information gathered above, we may reach the conclusion that the ReSL0 algorithm can reconstruct high-quality PAI images regardless of noise. Besides, the SL0 algorithm that does not include noise regularization terms cannot be used in a strong noise environment. The inequality optimization constraint in Eq. (11) can provide the better performance in the noise environment.

In Vitro Experiments.
Besides the numerical phantom experiments, the in vitro experiments were used to demonstrate the ReSL0 algorithm's performance. And the reconstruction results of the TwIST, SPGL1, and SL0 are also presented. A schematic diagram of the imaging system setup is displayed in Figure 7(a), where an experimental coordinate system ½x, y, z is also described. A Q-switched Nd:YAG pumped at 532 nm was used as the light source. The laser frequency is 10 Hz, and the pulse width is 6-10 ns. A focused piezoelectric transducer (Panametrics V309) with a central frequency of 5 MHz was controlled by a precision stepper motor for photoacoustic signal acquisition. The rotation radius of the transducer is 40 mm, and the rotation step size is 2°. At each sampling location, the photoacoustic signals were first amplified by a signal amplifier, then captured and averaged 30 times by a Tektronix MSO4000B mixed signal oscilloscope, and finally transmitted to the computer for signal processing and imaging.
The imaging sample used in the experiment is a gelatin cylinder with one graphite rod and two hairs as optical absorbers. Figure 7(b) is the photograph of the sample. The radius of the gelatin cylinder is 13 mm. The diameter of graphite rod is 0.5 mm, and the length of hair is 4 mm. In the phantom experiment, 60-view data and 90-view data are selected for reconstruction. The images are constructed by the TwIST, SPGL1, SL0, and ReSL0 algorithms, respectively.
The in vitro experiment results are displayed in Figure 8. The first and the second row show the reconstruction results from 60-view and 90-view experiment data. As can be seen from the first row of Figure 8, there will be a lot of noise in the reconstructed image and some low-contrast features will disappear when using 60-view experiment data. When the sampling data is sufficient, such as 90-view data, all four algorithms are able to construct good quality images. And the sizes and locations of the optical absorbers are well reconstructed. However, the resolution of reconstructed images using TwIST and SL0 algorithms is not as good as SPGL1 and ReSL0 algorithms. By comparing the results of numerical experiments and in vitro experiments, we can conclude that the quality of reconstructed images in vitro experiment is not as good as that of numerical experiments. Therefore, we need to collect more data to achieve exact reconstruction.

Conclusion
In this paper, we carried out and estimated a L 0 norm-based ReSL0 algorithm for the PAI. The main motivation of it is to replace the equality constraint with inequality constraint that allows some errors and make use of a regularization parameter to achieve a balance between the sparsity and the residual of objective function. The effectiveness and universality of this algorithm are demonstrated through numerical and in vitro experiments. Both visual inspection and quantitative measure comparisons have manifested that the ReSL0 algorithm can reconstruct better images than L1 norm and TV norm-based CS algorithms. Furthermore, the ReSL0 algorithm has similar computational efficiency to the SL0 algorithm and at the same time has better immunity to noise. Finally, the ReSL0 algorithm can significantly reduce the number of ultrasonic sensors and scanning time required to reconstruct high-quality PAI images.

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.