Photoacoustic image improvement based on a combination of sparse coding and filtering

Abstract. Significance: Photoacoustic imaging (PAI) has been greatly developed in a broad range of diagnostic applications. The efficiency of light to sound conversion in PAI is limited by the ubiquitous noise arising from the tissue background, leading to a low signal-to-noise ratio (SNR), and thus a poor quality of images. Frame averaging has been widely used to reduce the noise; however, it compromises the temporal resolution of PAI. Aim: We propose an approach for photoacoustic (PA) signal denoising based on a combination of low-pass filtering and sparse coding (LPFSC). Approach: LPFSC method is based on the fact that PA signal can be modeled as the sum of low frequency and sparse components, which allows for the reduction of noise levels using a hybrid alternating direction method of multipliers in an optimization process. Results: LPFSC method was evaluated using in-silico and experimental phantoms. The results show a 26% improvement in the peak SNR of PA signal compared to the averaging method for in-silico data. On average, LPFSC method offers a 63% improvement in the image contrast-to-noise ratio and a 33% improvement in the structural similarity index compared to the averaging method for objects located at three different depths, ranging from 10 to 20 mm, in a porcine tissue phantom. Conclusions: The proposed method is an effective tool for PA signal denoising, whereas it ultimately improves the quality of reconstructed images, especially at higher depths, without limiting the image acquisition speed.


Introduction
proposed method can be defined as a combination of low-pass filtering and sparse coding, we call it LPFSC.
The rest of this paper is organized as follows. Section 2 discusses the theory of TV denoising and ADMM methods, which are used in this paper. Section 3 describes the proposed denoising method, validation studies, and evaluation criteria. The experiments conducted to evaluate the performance of the proposed approach are described in Sec. 3, and the results are discussed in Sec. 4, followed by conclusions in Sec. 5.
2 Theoretical Background

Principles of Total Variation Denoising
Recently, the TV approach that promotes the sparsity of signals in the gradient domain has attracted significant attention in signal denoising applications. 48 The goal of TV denoising technique is to efficiently estimate and recover the desired N-point signal S ¼ fsðnÞg N n¼1 with the sparse or sparse-derivative representation from the measured noisy signal x, which is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 5 3 9 xðnÞ ¼ sðnÞ þ wðnÞ; (1) where wðnÞ is considered as additive Gaussian noise with the variance of σ 2 . TV denoising could be defined as the constrained minimization problem of a non-differentiable cost function in terms of the l 1 norm as below: where D as the first-order difference matrix is of size N × ðN − 1Þ and Ds is the first-order difference of an N-point signal sðnÞ. 49 With proper regularization parameter selection, Eq. (2) could be converted to the unconstraint problem as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 3 8 6 Tvdðx; λÞ ¼ arg min In this optimization problem, the regularization parameter λ plays a significant role and controls the degree of smoothing, so that when the λ ¼ 0, there is no smoothing and the result is the same as minimizing the sum of squares. On the other hand, increasing λ assigns a higher weight to the second term of Tvdðx; λÞ, which measures the oscillation of the desired signal sðnÞ and makes the solution sðnÞ piecewise. Although there are different algorithms to solve the TV denoising problem, majorization-minimization (MM) is found suitable to solve the optimization problem, which is hard to solve directly. 50

Alternating Direction Method of Multipliers
ADMM is a simple but powerful algorithm to solve a convex optimization problem by breaking it into smaller subproblems. 51,52 The ADMM algorithm is designed to solve the separable convex problems of the form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 1 9 0 min fðxÞ þ gðyÞ; subject to Ax þ By ¼ c; (4) where x ∈ R n , y ∈ R m , A ∈ R p×n , and B ∈ R p×m . The augmentation Lagrangian in Eq. (4) can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 4 0 L ρ ðx; y; λÞ ¼ fðxÞ þ gðyÞ þ λ T ðAx þ By − cÞ þ ρ where ρ is the penalty parameter, which is considered positive, and λ is the Lagrangian multiplier. Equation (5) is solved using three steps: x-minimization and y-minimization, which are split into N separate problems, and an updating step for multiplier λ as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 7 3 5 x kþ1 ≔ arg min x L ρ ðx; y k ; λ k Þ; y kþ1 ≔ arg min y L ρ ðx kþ1 ; y; λ k Þ; λ kþ1 ≔ λ k þ ρðAx kþ1 þ By kþ1 − cÞ: (6) 3 Method and Materials

Low-Pass Filtering and Sparse Coding
To solve the problem of PA signal denoising, the PA signal is modeled as the measured noisy signal xðnÞ, which is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 5 8 9 xðnÞ ¼ s lowfrq ðnÞ þ s sparse ðnÞ þ wðnÞ; where the desired signal sðnÞ includes two main components, a low-frequency component and a sparse or a sparse-derivative component. Here, s lowfrq ðnÞ represents the low-frequency component of the desired signal, s sparse ðnÞ represents the sparse or the sparse-derivative components of the desired signal, and wðnÞ is considered as additive Gaussian noise with the variance of σ 2 .
Since we are looking for efficient estimation of s sparse ðnÞ and s lowfrq ðnÞ, considering Eq. (7), in the following of s sparse ðnÞ estimation, one can estimate s lowfrq ðnÞ as follows: Considering the assumption that the frequency response of the low-pass filter is approximately zero-phase, we can conclude that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 3 3 3 highpass½xðnÞ −ŝ sparse ðnÞ ≈ wðnÞ; (11) where high pass refers to the high-pass filter. To achieve a computationally efficient approach, a zero-phase non-causal recursive high-pass filter which is proposed in Ref. 49 was deployed in our method. Since the PA signal is sparse and has a sparse derivative, the cost function of the optimization problem contains a linear combination of two regularization parameters, which promote piecewise smooth and sparse solutions. Therefore, the denoising problem, shown in Eq. (3), can be expressed as the unconstraint minimization problem of a non-differentiable cost function in terms of l 1 norm as below: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 2 1 7 arg min There are many solutions, such as MM and ADMM, to solve Eq. (12). Since the denoising process and image reconstruction speed are vital to achieve real-time PA imaging, we proposed to use the hybrid consensus ADMM method 53 to achieve a linear convergence and accelerating the conventional ADMM.

In-silico study
A simulation study generated by the k-wave toolbox in MATLAB ® (Mathworks, Massachusetts) was performed to evaluate the performance of LPFSC in PA signal denoising. 54 The initial pressure distribution is given by a 512 × 512 pixel image representative of a vascular structure. A 10-mm square grid was created, and the circular array detectors with 9-mm diameter and 60 elements, evenly spaced, were located around the region of interest (ROI) to receive the propagated PA wave from the object. We considered the center of each sensor as a point source. The sound speed was considered to be 1500 m∕s. The corresponding time array has 1019 data points that are 9.259 ns apart from each other (108-MHz sampling frequency). The input size of the PA signal was assigned 1019 × 60, and the reconstructed images have 128 × 128 pixels.

Experimental PA data acquisition setup
To further evaluation of the proposed PA signal denoising method and its effects on the quality of reconstructed PA images, phantom experiments were performed. The designed phantoms and the imaging setups are shown in Fig. 1. The first phantom contains two light-absorbing filaments with a diameter of 150 μm that were placed 1 mm apart from each other inside a water tank [ Fig. 2(b)]. An Nd: YAG/OPO nanosecond pulsed laser (Phocus core system, OPOTEK Inc., Carlsbad, California) with the pulse repetition rate of 10 Hz at wavelengths of 680 nm was used to illuminate the phantom. An ultrasound scanner (US) (Vantage 128TM, Fig. 1 Schematic of (a) the experimental setup used for the PA imaging of (b) two light-absorbing filaments that were placed inside a water tank and (c) cross-section view of blood-filled tubes embedded in a porcine tissue phantom.
The second phantom contains three polytetrafluoroethylene tubes (1-mm diameter), filled with human blood and embedded within a 30-mm-thick porcine tissue background. Bloodfilled tubes were placed inside the porcine tissue background at different depths from about 10 to 20 mm with 5-mm increments. PA acquisition was performed with the laser energy of 3 mJ∕pulse at the wavelength of 680 nm. A fiber bundle with a diameter of 20 mm used for guiding the laser light to the tissue. A 64-element phased-array US endoscopic transducer, with an active aperture of 9-mm long, was coupled to the phantom to acquire acoustic signals and provide high-resolution sector images [ Fig. 1(c)].

PA Signal Evaluation and Analysis
To evaluate the proposed approach, evaluation criteria such as peak signal-to-noise (PSNR) for signals, structural similarity index (SSIM), and contrast-to-noise ratio (CNR) for estimated images were used.
The PSNR as a common criterion to measure the quality of signal denoising based on the maximum possible value in the signal and mean square differences between denoised and reference signals is expressed in term of the logarithmic decibel scale (dB) as below: E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 3 . 3 ; 1 1 6 ; 6 4 4 PSNR ¼ 20 log 10 S max ffiffiffiffiffiffiffiffiffiffi MSE p where MSE is defined as a mean-square-error, S original and S denoised are original and denoised signals in size of M × 1, respectively, and S max is maximum possible value in signals.
We created PA images that represent an optical absorption distribution map of the targets via the conventional delay-and-sum (DAS) approach as the most commonly used reconstruction method in the PAI area. 55 For reconstructed images, the SSIM (in a scale of 0 to 1) as one of the most common criteria for image quality assessment and for evaluating the similarity of images (i.e., reference image and reconstructed images) is defined 56 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 3 . 3 ; 1 1 6 ; 5 1 0 where μ original and μ estimated are the mean of the original and estimated images, respectively, and also σ 2 original and σ 2 estimated are the variances of the original and estimated images, respectively. It is worth to mention that σ original;estimated is the covariance between the original and estimated image. The values of c 1 and c 2 are considered as constant values to avoid instability when the sum square of means or variances are very close to zero.
Finally, the CNR is a common criterion to determine image quality, especially in denoising processes with below definition: where S i and S O are the average intensity inside and outside of the objects, respectively. The σ O represents the standard deviation of the background. The background was defined as the pixels located inside the green dashed rectangular region selected in each set of PA images. For phantom studies, we considered the average of all frames as a ground truth (reference image) in each experiment.

Results and Discussion
To assess the proposed PA signal denoising method, we validated our method on numerical vessel phantom and experimental data of phantoms. The simulated PA signal, which was generated by the k-wave is shown in Fig. 2(a). Four different levels of additive Gaussian white noise with SNR levels 10, 5, −5, and −10 dB were added to the original simulated clean PA signal, and results of the proposed denoising method for the simulated noisy signals are shown in Figs. 2(b)-2(e). When the noise level was increased (for example, in SNR level −10), the PA signal peak was almost buried in added noise. As we have shown in Fig. 2, the LPFSC could suppress the noise and reconstruct the peak of the original signal, indicating the ability of the TV approach to recover the sparse or sparse-derivative signals. In the worst case (SNR level −10 dB), the difference between the original signal peak and the denoised signal peak is about 4%. This variance considering increasing the peak of noisy signal about 90% in comparison with the original signal is negligible. TV-based denoising is the most appropriate method for piecewise constant signals and preserves sharp edges in the underlying signal without requiring any step-size parameter as the amount of peak for the denoised signal. 57 It is worth mentioning that the proposed method could recover denoised signal until the SNR decreased to −15 dB. For further evaluation, the LPFSC signal denoising method was compared with two wellknown and widely used approaches of averaging and wavelet-based denoising methods. As shown in Fig. 3, the PSNR of the LPFSC method was compared to the frame averaging (using 20 frames) and wavelet-based denoising methods, at different noise levels. For PA signal denoising by wavelet method, a commonly used Symlet 6 wavelet with six-level of decomposition and Stein's Unbiased Risk Estimation threshold (SURE threshold) were selected. 29 The quantitative results obtained with the simulations show that the LPFSC compensates the low SNR of PA signal and outperforms the competing wavelet and averaging denoising methods in terms of PSNR by 24% and 26%, respectively, across all simulated noise levels.
The reconstructed images of denoised simulated PA signals with considering SNR of −10 dB for three methods of averaging, wavelet-based denoising, and LPFSC are shown in Fig. 4. The wavelet method has not been successful in recovering the small size vessels, which are corrupted by the noise. In addition, the averaging method using 20 frames cannot fully clean the signal due to presences of the coherent noise; however, our method using one frame could fully recover the original image and its details.
For the quantitative evaluation of the reconstructed images of denoised simulated PA signals, two SSIM and CNR criteria were used. Comparison between the results of three different methods, including wavelet-based signal denoising, averaging, and LPFSC in terms of SSIM and CNR with different levels of noise are shown in Table 1. On average (for different noise levels), the proposed denoising approach offers better CNR up to about 28% and 30%, and higher SSIM of about 23% and 24% in comparison with averaging and wavelet denoising methods, respectively. With increasing SNR levels from 5 to 20 dB and the reduction of noise level, the performance of all methods are improved. Although, for the lower SNR such as 5 dB, LPFSC outperforms wavelet and averaging denoising method by improvement about 58% and in terms of SSIM.
During experimental studies, we evaluated our method in a set of phantom studies in which two small-sized absorbers were placed inside a water tank, illuminated from the side and images from the top using a US transducer operating at frequencies between 4 and 11 MHz. The PA signal denoising results for the acquired signals of different detectors from the 20th frames of phantom data are shown in Fig. 5. The experimental results of PA signal denoising show significant improvement for the PSNR of PA signal of wires phantom about 35% by recovering peaks of original signal and reduction of noise.
Additionally, reconstructed images of denoised PA signals through three different methods are shown in Fig. 5 and are compared in Table 2, in terms of SSIM, CNR, axial and lateral full width half maximum (FWHM) for experimental phantom data. The selection of FWHM is to demonstrate the geometrical accuracy of reconstructed PA images. The experimental results of PA signal denoising prove that the LPFSC beats the performance of averaging denoising method (using all 20 frames) as the reference in the term of SSIM and CNR criteria by the 17% and 38%. The LPFSC method provides a closer to real size reconstructed PA image with a mean of axial FWHM of 0.77 AE 0.12 mm and lateral FWHM of 1.05 AE 0.04 mm, for two objects in experimental data. In comparison with the reconstructed image of one frame without any signal denoising with axial FWHM of 0.85 AE 0.18 and lateral FWHM of 1.08 AE 0.05 mm, our method could preserve lateral FWHM and improve axial FWHM about 10%. Whereas the wavelet denoising method did not improve the lateral and axial FWHM, and the averaging method was not effective for axial FWHM. It is worth to mention that both of these methods improve CNR of the reconstructed image in comparison with the reconstructed image of one frame with CNR 19.32 dB. Therefore, one of the main improvements gained by the proposed LPFSC method is having reduced geometrical distortion as well as a high contrast at the same time. More importantly, the LPFSC method offers using one frame in comparison with the averaging method, which requires 20 frames. Using the hybrid ADMM, which stands out as efficient and easily implementable, leads to fast convergence of our method.
Finally, we evaluated the performance of the proposed denoising method with data acquired by an endoscopic probe and from the porcine tissue phantom. Since there was no significant  difference between the results of frame averaging and wavelet method in in-silico study, the reconstructed PA images of denoised signal with frame averaging and LPFSC methods are compared in Fig. 6. The images form by averaging of 10 frames were used as a ground truth or the reference image. All constructed images from denoised signal using LPFSC and traditional frame averaging with the same number of frames (n ¼ 1 to 7) were compared to the ground-truth reference to calculate the SSIM. The number of averaged frames varied between one and seven. First and second rows related to the results of the averaging method and our proposed method, respectively. The walls of blood-filled tubes are generating strong PA signals that are visually presented as two parallel bright lines in PA images. To calculate CNR, an ROI located within the tube [as shown in Fig. 6(e)] was used to measure the signal from the target. To further quantify the image enhancements, CNR and SSIM of frame averaging and LPFSC using five frames are compared in Fig. 6(h). The use of five frames averaging in conjunction with LPFSC markedly reduced the background noise, which improved the mean of CNR and SSIM by 32% and 31%, respectively.  Evaluation of CNR and SSIM for different numbers of averaged frames and at different imaging depths are shown in Fig. 7. As anticipated, CNR and SSIM parameters are improved with increasing the number of averaged frames. However, the proposed denoising method provides higher CNR and SSIM using the same number of averaged frames for objects placed at different depths. The mean of CNR using LPFSC with one frame was improved 63% in comparison with averaging for the objects evaluated at three different depths. A closer look at the results reveal that the mean CNR using the proposed denoising method using three averaged frames is higher than averaging only and using seven averaged frames. This comparison clearly shows the potential of using the proposed method to enhance the quality of PA images with a smaller number of averaged frames, which can lead to a preserving image quality at higher imaging speed. Our results indicated the LPFSC with three frames averaging is capable of detecting the objects located at depth 3 with SSIM of 0.91, which shows 90% improvement compared to the standard averaging method. The average measured SSIM for objects located at three different depths was calculated as 0.88 AE 0.15 when LPFSC is utilized. This result indicates 33% improvement compared to averaging alone.
In this study, a range of cutoff frequencies was experimentally tested on data sets. The accepted range of cutoff frequencies was concluded to be between 0 and 0.3 cycles sample for different data sets, experimentally. For in-silico and first experimental phantom studies, the optimal value of cutoff frequency was selected 0.1 cycles sample and for porcine tissue data, it was selected 0.15 cycles sample . Also, we consider equal weight for the TV and LPF; therefore, the lambda value was selected equal to be 1, experimentally. Changing lambda within the range of 0.8 to 1.2 showed acceptable results for PA signal denoising. For lambda smaller than 0.8, the amplitude of denoised signal is decreased compared to the original signal. The lambda bigger than 1.2 downgraded the performance of LPFSC. The performance of LPFSC could also be affected by the adjustment of low-pass filter cutoff frequency. However, small deviations from the ideal cutoff frequency up to 20% could be compensated by the sparse denoising part of the proposed algorithm.

Conclusion
The efficiency of PA imaging is routinely limited by the presence of background noise and suffering from low SNR, which resulted in the poor quality of the reconstructed images. Since the PA signals can be modeled as a sum of two low frequency and sparse components, we proposed a denoising approach that simultaneously estimates a low-pass and a sparse signal from an acquired noisy signal based on TV optimization approach, and using hybrid ADMM. Both in-silico and experimental work on tissue mimics were used to evaluate the performance of the proposed technique. The results demonstrated that the LPFSC method possess a superior performance to compensate the low SNR PA signals and offered a better CNR and SSIM for the reconstructed images compared to the frame averaging method. This comparison clearly shows the potential of using our proposed method to enhance the quality of PA images while maintaining high-speed imaging which is an essential need in many applications of PA imaging. In other words, the frame rate of PA images, which is always a challenge in real-time PA imaging can be significantly improved.

Disclosures
The authors declare that they have no conflicts of interest.