1. Introduction

Photoacoustic tomography (PAT) is a hybrid imaging modality providing optical absorption contrast at high ultrasonic resolution [1–4]. It involves usage of a pulsed laser (temporal duration in nanoseconds) in the range of 600 - 1000 nm (near-infrared (NIR)) for irradiating the tissue. The absorbed light by the tissue chromophores tends to increase the temperature which leads to thermoelastic expansion. The photoacoustic (PA) waves are generated because of light absorption and propagate inside the tissue, to be acquired by the transducers placed at the boundary of tissue. The captured measurements by transducers are then used as an input to the various reconstruction algorithms for estimating the initial pressure rise distribution inside the tissue. The strength of PAT lies in the fact that the attenuation and scattering of acoustic waves is far less compared to light propagation and thus the propagation inside the tissue without significant scattering can be achieved. Photoacoustic imaging has been widely used for revealing functional and structural information in clinical and pre-clinical applications, non-invasive monitoring of traumatic brain injury [5], oncology [6,7], pathology, molecular imaging [3] and enables deep tissue monitoring because of higher light penetration in the NIR-window [8].

One of the bottlenecks for translation of PAT to pre-clinical/clinical applications is the lack of robust reconstruction methods being available for generating accurate images in photoacoustic tomography. The present analytical methods for computing the result of the inverse problem are composed of filtered backprojection and Fourier transform based methods [9,10]. Analytical algorithms are generally computationally efficient but require large number of ultrasonic sensors around the imaging domain. These algorithms also require more number of ultrasonic transducers or acquiring more number of time samples, which increases the data acquisition time. In general, the analytical or time reversal based methods have limited quantitative accuracy [11–13].

The term quantitative in here refers to the reconstructed initial pressure distribution being more close to the true pressure distribution in terms of contrast recovery. Since the reconstruction step is performed with limited number of transducer elements and limited number of time samples, it results in inferior contrast recovery and in here we deal only with the acoustic inverse problem. The model based techniques have gained wide importance in the recent past to solve the inverse problem and get quantitatively accurate images for limited data cases [9,14]. The term limited data refers to the limited amount of the measurement data being acquired utilizing photoacoustic experimental setups. There are clinical experimental setups with the number of acoustic sensors being as high as 512 and the number of time samples could be 2048 [15]. In this work, the number of time samples used were only 512 and the maximum number of acoustic detectors being only 100, thus being limited data in nature. The limited data case utilizes comparatively less number of detectors than the full data cases. Here, we have utilized 100 detectors as an example, but the detectors can be decreased further to show the utility of the proposed method. Here, we have not included the limited view cases which cover only part of the simulation grid as the focus is only on limited data cases.

The conversion efficiency of light to sound in photoacoustic effect is quite low and the signal to noise ratio (SNR) of photoacoustic signals in case of thick tissues is also limited, many methods have been proposed to improve the SNR of these signals. In Ref. [16], a method was proposed combining the empirical mode decomposition and mutual information to reduce noise in photoacoustic imaging. Another work introduced exponential filtering of the SVD coefficients to improve the reconstruction and mitigating the effects of the modeling errors present during the reconstruction [17]. Even application of total least squares to compensate the modelling errors and tackle the low SNR of raw data was also investigated [18]. In Ref. [19], a SVD based model was introduced for analysis of photoacoustic imaging system and utilized denoised imaging operator for estimating the number of measurable singular vectors for the system. Another investigation to identify and remove laser noise in photoacoustic images utilizing an ultrasound clinical scanner using SVD was also carried out [20]. Compressed sensing (CS) and deep learning (DL) based methods for cases of sparse data were also applied to the problem at hand [21–23]. Various other methods, such as $l_1$ minimization, total variation have also been utilized as a sparsity regularizer for denoising as well as reconstruction scenarios [24]. Total Variation (TV) measures the gradient of the image for getting the variations in an image and has been shown as an effective regularization method for solving linear inverse problem [25]. In Ref. [26], an adaptive steepest descent projection onto convex sets (ASD-POCS) involving the TV step for limited data reconstruction was investigated. In similar line, Ref. [27] proposed a gradient descent algorithm based on TV for undersampled measurements, and Ref. [28] proposed a TV regularization enhanced by bregman iterations for solving the PAT sub-sampling problem for increasing the acquisition speed and achieving good resolution. In Ref. [29] deep learning was applied to improve the reconstruction result obtained by truncated SVD. On the data side, a deep learning based method was proposed for super-resolving the sinogram (photoacoustic raw data) from undersampled measurements to obtain super-resolved, denoised and bandwidth enhanced sinogram for improving the photoacoustic image reconstruction [30]. In Ref. [31], a deep learning method based on learned regularizer was proposed for mitigating the aliasing artifacts obtained in limited data photoacoustic problem.

The model based regularization techniques are more robust for noisy cases and the reconstruction obtained is quantitative in nature. Several model based methods, including Tikhonov based regularization, promote smoothness in the resultant reconstructed photoacoustic image. In simple words, the result is generally blurred and advanced computational methods utilize deblurring techniques, such as Basis Pursuit Deconvolution (BPD) [32] to provide sharper images. All these methods result in a solution which is biased to the choice of regularization parameter. Regularization methods such as total variation based regularization which minimizes the variation in the solution can provide sharper features in the reconstructed image, but also requires careful choice of regularization parameter [33]. Generally noise is present in the measurement data which necessitates regularization to mitigate the effect of noise and often this regularization can induce smooth features (blurriness) in the reconstructed photoacoustic image. The iterative methods are effective in removing the noise as compared to the analytical techniques and getting a reconstruction close to the ground truth and hence provide more meaningful solutions.

Recently proposed method, framed in Plug and Play (P & P) framework, uses the existing denoising algorithms for getting a solution of various tasks particularly to improve and accelerate the image reconstruction [34]. These P & P methods propose a way to decouple the image prior with the measurement model. The image prior is handled using a standard denoising operation such as total variation etc. The major advantage of using such a model is that the prior is not defined explicitly (it is implicitly defined by the denoiser [35,36]). The technique of P & P has found many application areas as in post-processing of compressed image [37], bright field electron tomography [38], and Poisson denoising [39]. In Ref. [40] another technique was proposed to reduce the number of iterations as well as less parameter tuning to achieve the same. The cost function was transformed into a novel optimization problem and solved using efficient techniques to get the solution. The problem with this technique is that for large matrices the inverse needs to be calculated explicitly, which is a time consuming process having high computational complexity. Also it is not always possible to get an inverse for an ill-conditioned matrix. The computation of a pseudo-inverse often requires the use of regularization parameter e.g. in case of SVD it is the number of singular values to get the accurate pseudo-inverse. Here a SVD [41] based form of the Iterative Denoising and Backward Projection based method was utilized to get the pseudo-inverse.

The proposed Singular Value Decomposition based Iterative Denoising and Backward Projection (SVD-IDBP) method is a novel technique to include any state of the art denoising technique to be utilized as priors for solving the inverse problems and can be easily included in all inversion schemes without any overhead [34]. Here, total variation denoising was deployed as a denoising prior, but other models such as Non local Means [42], Guided filter [43], and wavelet thresholding (soft and hard) can also be seamlessly utilized. The proposed technique is agile and effortless to be integrated into the reconstruction framework (software’s) to provide improved solution of inverse problems, examples include image reconstruction and image inpainting [34].

In the proposed Plug & Play based method, the optimization problem is split into two parts - one is the prior term where the denoising is involved using state of the art denoisers and the other is the forward model used for getting the sinogram domain data. These are decoupled (solved independently) and hence provide improved solution to the inverse problem. The benefit of decoupling is that one branch composed of denoising is directly dependent on the prior, while the other branch only depends on the forward model comprising the model based reconstruction using SVD and $l_2$ regularization. The derived SVD based inversion in the model also helps in computation of the pseudo-inverse for ill-conditioned matrices (common in ill-posed inverse problems like photoacoustic tomographic imaging).

Specifically, the SVD [41] based reduced dimension plug and play model for improving the image reconstruction for photoacoustic tomography was introduced in this work. This iterative SVD based model can be used for solving ill-conditioned linear equations, where the direct computation of inverse is not plausible [44]. The advantages of using the SVD based iterative plug and play model (SVD-IDBP) are as follows:

2. Photoacoustic (PA) image reconstruction

The acoustic waves that are generated due to heating of the tissue by irradiation of laser source are modeled by the photoacoustic wave equation. The forward model computes the pressure waves which are captured by the transducers placed at the boundary. These PA waves that propagate inside the tissue follow the equation: [4]

Here, $P(z,t)$ : the pressure at position $z$ and time $t$, $c$ : the speed of sound, $C_p$ : specific heat capacity, $\beta$ : the thermal expansion coefficient, and $E(z,t)$ : the energy deposited per unit volume per unit time.

2.1 Building the system matrix

The numerical framework utilized in here for the computation for forward and backward problems here is similar to previously used one (available in Refs. [ [46]]). Here it is reviewed briefly for completeness. The dimensions of the imaging domain utilized here for the framework is ${n \times n}$ which is reshaped to make a long vector of dimensions ${n^2 \times 1}$. The dimensions of the obtained system matrix $\textbf {A}$ is ${m \times n^2}$. Here the variable ${m}$ represents the product of the total number of ultrasonic detectors and the number of time samples acquired. The response of each pixel is computed using the forward model and is stacked as a column vector to make the complete system matrix. The measurement vector obtained by collecting the data using the detectors at the boundary is denoted as ${y}$ having dimension of ${m \times 1}$ [46]. The forward model for the photoacoustic imaging is given as

2.2 Linear back projection (LBP) Method

The simplest linear back projection (LBP) output can be obtained as [32,47]

Here $\textbf {T}$ denotes the transpose of the matrix. The backprojection technique is only qualitative in nature and hence not used for quantitative imaging [14]. The backprojection based output is generally used as the initial point for starting the iterations for the iterative methods [46].

2.3 Lanczos Tikhonov regularization method

The inverse problem arising in photoacoustic tomography is ill-conditioned for limited data (the number of sensors used are comparatively less) case. In general, Tikhonov regularization has been popular for solving this ill-posed problem of photoacoustic imaging. The plot of the singular values of the system matrix is shown in Fig. 1. A large fraction of the obtained singular values are close to zero owing to the ill-condition of the matrix. The solution obtained is computationally expensive as it takes ${\mathcal {O}}(n^6)$ operations. To reduce the computational complexity of this method, Lanczos bidiagonalization for the above problem was proposed [48]. This method has been previously used for providing optimal solution and has been the standard method to be deployed for solving the inverse problem of photoacoustic tomography [18,32,48]. For more details, the readers are encouraged to refer to Ref. [48].

 

Fig. 1. Example plot of the magnitude of singular values of the corresponding system matrix.

Download Full Size | PDF

2.4 Basis pursuit deconvolution (BPD) method

Basis Pursuit Deconvolution (BPD) [32] was employed as an additional post-processing method to deconvolve and get an improved solution as compared to the one obtained using the Lanczos Tikhonov regularization method. The regularization often blurs the obtained pressure distribution and thus the resulting regularization effect can be minimized using the BPD method. It deblurs the solution $x_{LTH}$ obtained using the Lanczos Tikhonov Regularization method. This is achieved by utilization of Split Augmented Lagrangian Shrinkage Algorithm (SALSA) [36] algorithm for minimizing of the cost function, which typically uses $\ell _1$-norm based regularization for promoting sharp features in the reconstructed solution (when the TV term in Algorithm 1 replaced with $\ell _1$-norm, it will lead to the form utilized in Ref. [32]).

2.5 Total variation regularization method

The Tikhonov regularization promotes smooth features in the reconstructed image due to the presence of the $l_2$ norm of regularization part. Other standard method for solving this ill-conditioned problem is to minimize the variation obtained in the solution along with the residual term. This method comprising the variation and the residual term is called as total variation (TV) regularization technique [49,50]. The function to be minimized for the total variation case can be written as

Here, ${\eta }$ denotes the regularization parameter and ${\|.\|_{TV}}$ represents the isotropic total variation [51].

Here, $\phi$ denotes the TV functional ${\|.\|_{TV}}$ and Alternating Direction Method of Multipliers (ADMM) framework was utilized here for TV regularization which has the same form as given in Ref. [52]. The Split Augmented Lagrangian shrinkage Algorithm (SALSA) (given in Algorithm-1), was utilized in this work. The LBP reconstruction was utilized as an initial guess of ${x}$ and the variables $l_0$ and $w_0$ were initialized to zero for the SALSA algorithm.

2.6 Singular value decomposition based iterative denoising and backward projection (SVD-IDBP) method [proposed]

The forward problem of the photoacoustic imaging is given in Eq. (2). If noise is present in the data, the equation is modified to

Here ${\|.\|_2}$ represents the ${l_2}$ norm of the vector and ${\tilde {x}}$ is the solution vector. The cost function is formulated in with a small change for deriving the proposed algorithm:

Knowing ${\| x \|^2_{ \textbf{A}^{\textbf{T}}\textbf{A}} \overset {\Delta }{=} x^T \textbf{A}^{\textbf{T}}\textbf{A}x}$ and $\textbf {(AB)}^\textbf{T}=\textbf{B}^\textbf{T} \textbf{A}^\textbf{T}$, the term $\|(\textbf{A}^{\mathbf{\dagger}}y-\tilde {x})\|_{\textbf{A}^{\textbf{T}}\textbf{A}}^2$ can be rewritten using these properties as

Using the SVD decomposition [41] of matrix $\textbf {A}$, the matrix $\textbf {A}$ can be decomposed as

The pseudo-inverse of $\textbf {A}$ can be written as:

Using Eqs. (12) and (13), the Eq. (14) can be written as

Thus, the product of pseudo-inverse and the matrix can be simplified using Eq. (12) as

Thus Eq. (7) can be modified as

This cost function of minimizing ${f({\tilde {x}})}$ over ${\tilde {x}}$ can be also written as

Similar to Ref. [53], here in order to get ${\tilde {x}}$, freedom is required on how to chose ${\tilde {b}}$. Hence we relax the above term by relaxing the constraint ${\tilde {b}= \textbf {V S}^{-1} \textbf{U}^{\textbf T} y}$ to ${y = \textbf {USV}^\textbf{T} \tilde {b}}$ such that no inverse is involved in this term. Since ${\tilde {b}}$ is always multiplied by $\textbf {A}$, the null space of $\textbf {A}$ is always ignored. Thus they will disagree with the prior ${p(\tilde {x})}$ as they are not controlled by these terms. To solve this problem, the term ${\frac {1}{2\sigma _{e}^{2}} \| \tilde {b}-\tilde {x} \|_{ \textbf {VS}^{\textbf T}\textbf{S V}^{\textbf{T}}}^{2}}$ is replaced with another term ${\frac {1}{2(\sigma _{e}+\delta )^{2}} \| \tilde { {b}}-\tilde { {x}} \|_{2}^{2}}$, where ${\delta }$ is another regularization parameter. Higher value of ${\delta }$ tends to reduce the effect of the fidelity term, while lower value of ${(\delta +\sigma _e)^2}$ will penalize the solution more. The problem of choosing an optimal ${\delta }$ is discussed in Ref. [53]. Thus the Eq. (19) is modified to

This equation can be solved in an iterative manner by using Alternating Direction Method of Multipliers technique. The ${\tilde {x}_k}$ is obtained using

This equation tries to find the solution ${\tilde {x}_k}$ with a denoiser for the white Gaussian noise having variance ${\sigma ^2 = ({\sigma _e+\delta })^2}$. This denoiser is applied iteratively on the input image ${\tilde {b}_{k-1}}$ and can be written in the compact notation as given

This equation defines a projection of ${\tilde {x}_k}$ onto the subspace and has a closed form solution

Thus the solution of the Eq. (19) can be written as

Using this solution and plugging it in Eq. (19) will result

In the original equation (Eq. (6)), the fidelity term tries to provide better fit between the measurements ${y = \textbf {A}\tilde {\textbf {x}}+e}$ and $\textbf {A}\tilde {\textbf {x}}$. Here the new problem aims at better fit between the terms

Here $\textbf {P}_{\textbf{A}} \overset {\Delta }{=} \textbf{A}^{\dagger }\textbf{A}$ is the orthogonal projection onto the row space of $\textbf {A}$. Assuming that the matrix $\textbf {A}$ is of full low rank, the matrix $\textbf {P}_\textbf{A}$ and $\textbf{A}^{\textbf{T}}\textbf{A}$ have rank $m$, even though they may have very different eigen values. The eigen values of $\textbf {P}_\textbf{A}$ can be only 1 in the row space of $\textbf {A}$ while 0 in the null space of $\textbf {A}$.

It is a well known that if the singular values of the linear operator are not spread over a wide range, then the linear least square optimization methods perform better [54]. Thus if the prior ${p(\tilde { {x}})}$ gives a strong restriction on $(\textbf {I}_{\textbf {n}}-\textbf {P}_{\textbf {A}}){\tilde {x}}$ given $\textbf {P}_\textbf {A}{\tilde {x}}$, then solving ${c\|{\textbf{A}^{\mathbf{\dagger}}y-\textbf {P}_\textbf {A}\tilde {x}\|_2^2} + s(\tilde {x})}$ might be more stable than solving ${c{\|y-\textbf {A}\tilde {x} \|_2^2} + s(\tilde {x})}$. When the noise is low and good image priors are highly non convex, a numerical optimization process for

Thus, the replacement provides a solution which is more close to the actual solution. The proposed SVD based method eliminates the need of the prior function ${p(x)}$. The prior function is chosen depending on the denoiser ${D(.;\sigma )}$. Here total variation denoiser was used and combined with the SVD based method for improved performance. As $D$ is an operator involving the image with the threshold parameter, can be used as a black box for denoising, other popular techniques such as Non local Means [42], Guided filter [43], and wavelet thresholding (soft and hard) can also be easily deployed in this framework. The output ${\tilde {b}_{k}}$ is a better approximation to the ground truth signal as compared to the raw observations ${y}$ obtained. Thus it iteratively switches between reconstructing the signal and iteratively using this reconstruction to improve. Overall, it provides better estimates of the measurements. The proposed algorithm, which was termed as Singular Value Decomposition based Iterative Denoising and Backward Projections (SVD-IDBP), was presented in Algorithm 2. The benefits of the proposed algorithm are that the offset can be chosen for the singular values below which they are assumed to be zero (or noise floor). It helps in removing the noise which was assumed to be present in the high frequency components and in the lower singular values of the SVD matrix. For very large matrices the effect of these noisy singular values can be mitigated using the proposed method as seen in the noisy cases (20 dB SNR) presented in this work [55].

3. Figures of merit

For comparing the performances of the reconstructed initial pressure distributions in the numerical phantom cases obtained using the proposed technique and other standard techniques, Root Mean Square Error (RMSE), Pearson Correlation (PC) and Contrast to Noise Ratio (CNR) were utilized. Since ground truth was not available for the experimental phantom and in-vivo case, the criteria used for comparison is the Signal to Noise Ratio (SNR). These figures of merit are defined as given below.

3.1 Root mean square error (RMSE)

It is a commonly used metric to predict the accuracy of the proposed method. It is defined as [56,57]

3.2 Pearson correlation (PC)

PC (range -1 to 1) measures the correlation between the reconstructed image and the target image and is defined as follows [58,59]:

3.3 Contrast to noise ratio (CNR)

CNR is defined as follows [15,60–62]:

3.4 Signal to noise ratio (SNR)

SNR of the reconstructed image is defined as [63]:

4. Numerical and experimental studies

Figure 2 shows the schematic of the acquisition geometry used for the collection of photoacoustic data along with the hundred ultrasonic/acoustic transducers which were positioned at the border of the domain to be imaged. The computational grid comprises of 100 ultrasonic/acoustic detectors, which were positioned at a length of 22 mm from the center of the imaging domain. The computational grid used for the imaging had a dimension of 50.1 mm ${\times }$ 50.1 mm and the reconstruction grid had dimension of 20.1 mm ${\times }$ 20.1 mm. The computational grid is the complete grid comprising the data generation and the reconstruction grid including the ultrasound transducers. The photoacoustic data for these was generated on a grid having dimension of 401 ${\times }$ 401, while a grid of dimensions 201 ${\times }$ 201 was utilized for performing the reconstruction of the photoacoustic image data. The dimension of the grid was same for generation and reconstruction, but the grid lines are changed as 401 ${\times }$ 401 for data generation while 201 ${\times }$ 201 for data reconstruction for mimicking the real experimental scenarios. Here, only the ring-scanning setup was utilized and other geometries for data acquisition including handheld probes such as linear arrays, curvilinear arrays, or hemispherical arrays are not utilized in this work [64,65].

 

Fig. 2. Schematic diagram depicting geometry of photoacoustic data acquisition. The position of the hundred ultrasonic transducers are shown by dots around the imaging domain. The size of the computational imaging grid is 50.1 mm $\times$ 50.1 mm while the size of the data reconstruction grid is 20.1 mm $\times$ 20.1 mm.

Download Full Size | PDF

Two different numerical phantoms having different characteristics were utilized for comparison of the proposed reduced dimension plug and play method with other standard techniques. Figure 3(a) shows the blood vessel phantom used which consists of thick and thin blood vessels. A phantom having letters ‘PAT’ with sharp edges (Fig. 3(b)) was also utilized. The performance of the proposed method was investigated by utilizing photoacoustic data with varying SNR levels from 20 dB to 60 dB and random white Gaussian noise was added for obtaining the prescribed SNR. The levels of 20 dB, 40 dB and 60 dB SNR corresponds to approximately 10%, 1% and 0.1% of the amount of noise present in the measurement data. Noise level of 40 dB SNR (1% noise) is the normal noise level generally present in the measurement data. The experiments were also performed for higher noise levels of 20 dB SNR (10%) and improved results were obtained as compared to the other state of the art methods. Previous studies have shown that the noise levels of 7% in the measurement data [29,31] being taken as the limit. Since at 20 dB SNR, the noise level is already 10% (representing the worst-case scenario), the studies pertaining to the proposed method failing with increasing further noise level were not performed. The ultrasonic transducers positioned on the boundary have a center frequency of 2.25 MHz with 70 ${\%}$ bandwidth were utilized in collection of this raw data. The photoacoustic data was collected at a sampling rate of 20 MHz and the number of time samples acquired by the transducers were 512. The dimension of system matrix $\textbf {A}$ is 51200 ${\times }$ 40401. The assumption is that system is homogeneous, non absorptive and non-dispersive. The $k-$ wave toolbox [66] was utilized for the generation of data for all the numerical simulations.

 

Fig. 3. The target numerical phantoms that were utilized in this work consists of the following: (a) Blood Vessel phantom and (b) ‘PAT’ phantom.

Download Full Size | PDF

The sinogram of the in-vivo mouse brain and the horse-hair phantom (that has triangular shape) has been acquired using Nd:YAG laser setup. This photoacoustic data obtained using the experimental setups has been used for knowing the utility of the proposed method in real time scenario. The details for the experimental setup used for acquiring data of the triangular-shaped horse hair phantom and the in vivo data was provided in Ref. [67]. It is pertinent to mention that the experiments conducted on animals as part of this work, the guidelines and regulations accepted by the institutional Animal Care and Use committee of Nanyang Technological University, Singapore (Animal Protocol Number ARF-SBS/NIE-A0263) were followed. Here, the reconstruction grid has a size of 40 mm ${\times }$ 40 mm containing 200 $\times$ 200 pixels. The photoacoustic data was acquired at a sampling frequency of 25 MHz (1024 samples), and subsampled at half to get 512 samples.

5. Results and discussion

The initial pressure distribution using the LTH method for the blood vessel phantom is shown in Fig. 4(a) for 20 dB SNR in the data. The reconstruction using Basis Pursuit Deconvolution (BPD) and Total Variation are shown in Fig. 4(b) and (c) respectively for the same SNR in the data. The proposed method result is shown in Fig. 4(d). The figures of merit (RMSE, CNR, and PC) discussed in this work for these results are shown as a bar plot in the first row of Fig. 6. The proposed reconstruction method improved the RMSE, CNR and PC by 14.15 ${\%}$, 254.80 ${\%}$ and 112.79 ${\%}$ respectively for the 20 dB SNR case. For low noise level the improvement obtained is very high as can be seen from the Fig. 4. Similarly for the SNR of 40 dB, the proposed method gave an improvement of 26.65 ${\%}$, 14.15 ${\%}$ and 3.34 ${\%}$ in RMSE, CNR and PC as compared to the other reconstruction methods. The figures for the LTH, BPD, TV and the proposed method are given in Figs. 4(e), (f), (g) and (h) respectively for the 40 dB SNR of the data. The artifacts are significantly reduced for this case. For SNR of 60 dB, the improvement using proposed method is relatively lower and is 8.51 ${\%}$, 7.01 ${\%}$ and 0.85 ${\%}$ in RMSE, CNR and PC as compared to the other reconstruction methods. The corresponding reconstruction results are presented in the last row of Fig. 4.

 

Fig. 4. Compares various methods of photoacoustic image reconstruction detailed in this work when applied on the blood vessel numerical phantom (Fig. 3(a)) data. The comparison was preformed for three different SNR levels of sinogram mentioned against each row (20, 40, and 60 dB). The reconstructed photoacoustic image using (a) LTH based method (Sec. 2.3) (b) Basis Pursuit Denoising (BPD, Sec. 2.4) (c) Total Variation (TV) regularization (Sec. 2.5) (d) proposed Singular Value Decomposition based Iterative Denoising and Backward Projections (SVD-IDBP) method (Sec. 2.6) for the SNR of data being 20 dB. Similar results were shown for 40 dB and 60 dB SNR cases in the middle and last rows. The bar plots of the figures of merit for the obtained results were given in the first row of Fig. 6

Download Full Size | PDF

The performance of the proposed reconstruction method in case of ‘PAT’ phantom case is shown in Fig. 5 for varying SNR levels in the data. The figures of merit for these reconstruction results were presented in the last row of Fig. 6. The same trend as observed in previous two cases was recorded here as well. The proposed reconstruction improved the RMSE, CNR and PC by 67.98 ${\%}$, 523.42 ${\%}$ and 135.33 ${\%}$ respectively for the 20 dB SNR case. Similarly for the 40 dB SNR case, the proposed method gave an improvement of 15.36 ${\%}$, 11.89 ${\%}$ and 0.93 ${\%}$ in RMSE, CNR and PC as compared to the other reconstruction methods. For 60 dB SNR case, the improvement was 10.59 ${\%}$, 15.08 ${\%}$ and 1.30 ${\%}$ in RMSE, CNR and PC as compared to the other reconstruction methods.

 

Fig. 5. Compares various methods of photoacoustic image reconstruction detailed in this work when applied on the ‘PAT’ numerical phantom (Fig. 3(c)) data. The comparison was preformed for three different SNR levels of sinogram mentioned against each row (20, 40, and 60 dB). The reconstructed photoacoustic image using (a) LTH based method (Sec. 2.3) (b) Basis Pursuit Denoising (BPD, Sec. 2.4) (c) Total Variation (TV) regularization (Sec. 2.5) (d) proposed Singular Value Decomposition based Iterative Denoising and Backward Projections (SVD-IDBP) method (Sec. 2.6) for the SNR of data being 20 dB. Similar results were shown for 40 dB and 60 dB SNR cases in the middle and last rows. The bar plots of the figures of merit for the obtained results were given in the first row of Fig. 6

Download Full Size | PDF

 

Fig. 6. The bar plots for figures of merit, RMSE, CNR, and PC, are shown for the proposed method and the state of the art techniques: (a-c) for blood vessel phantom (Fig. 3(a)), (d-f) for ‘PAT’ phantom (Fig. 3(b)).

Download Full Size | PDF

Overall, in the numerical phantom cases, the improvement with utilization of plug and play priors was significant (as high as 5 times in terms of figures of merit) especially in cases of noise level being high (SNR of data being low). For the low noise case (SNR being close to 60 dB), the plug and play priors did not improve reconstructed image quality significantly and performed on par with TV regularization method.

Figure 7 shows the reconstruction results obtained in the case of experimental phantom of triangular shaped horse hair. Figures 7(a), (b) and (c) show the reconstruction result using the LTH based method, BPD and TV respectively. The reconstruction obtained using the proposed method is shown in Fig. 7(d). The corresponding figure of merit (SNR) was given below each of this result. The proposed reconstruction method improves the reconstruction by reducing the artifacts as well as improving the ${{SNR}_r}$ by 8.1 dB. Here, the main aim was to show the utility of the iterative SVD based plug and play priors for improving the photoacoustic image reconstruction. As an extension to this work in future, the same technique can be applied for phantoms having different absorption coefficients. Similarly, Figs. 8(a), (b) and (c) show the reconstruction using LTH based method, BPD and TV respectively for the in vivo mouse brain data. The reconstruction obtained using the proposed method is shown in Fig. 8(d). The improvement in terms of figure of merit ${{SNR}_r}$ is 4.21 dB in comparison to TV-regularization based method (which provides the best result among standard methods discussed in this work) for this case. The results obtained using any iterative method depend on the regularization parameter, which can be chosen using the GCV [68], L-curve [69] or another optimization technique [70,71]. These methods are computationally demanding [71]. As the main aim of the work was to show the utility of the iterative SVD based plug and play priors over other iterative techniques and hence the regularization parameter was chosen heuristically and tuned to give the best performance for any algorithm. The number of steps of bidiagonalization were chosen as 40 and $\alpha$ (regularization parameter) as 0.3, while the regularization parameter for the BPD method was chosen as 1e-5 (and the maximum no of iterations are kept as 10000) and are obtained in the same way as in Ref. [32]. The regularization parameters for the TV based regularization was 1e-3 and the rest parameters were chosen utilizing the technique proposed in Ref. [36,72]. The regularization parameters for the proposed SVD-IDBP method were chosen as $\beta$=0.5, lagrange parameter as 1, and regularization parameter as 0.018 for the TV denoising. The regularization parameters for the experimental reconstruction was also chosen in the same way as for the numerical phantoms. The parameters tuned for one numerical experiment case were utilized in rest experiments to avoid fine tuning of these parameters for every case.

 

Fig. 7. The reconstruction of initial pressure distribution in experimental horse hair phantom with utilization of (a) Lanczos Tikhonov Regularization (Sec. 2.3), (b) Basis Pursuit Denoising (BPD, Sec. 2.4 ),(c) Total Variation (TV) regularization (Sec. 2.5), and (d) proposed Singular Value Decomposition based Iterative Denoising and Backward Projections (SVD-IDBP, Sec. 2.6). The corresponding figure of merit, SNR of the reconstructed image, is provided below each reconstructed image.

Download Full Size | PDF

 

Fig. 8. The reconstruction of initial pressure distribution in experimental in-vivo mouse brain with utilization of (a) Lanczos Tikhonov Regularization (Sec. 2.3), (b) Basis Pursuit Denoising (BPD, Sec. 2.4 ),(c) Total Variation (TV) regularization (Sec. 2.5), and (d) proposed Singular Value Decomposition based Iterative Denoising and Backward Projections (SVD-IDBP, Sec. 2.6). The corresponding figure of merit, SNR of the reconstructed image, is provided below each reconstructed image.

Download Full Size | PDF

Typical computational time recorded for the methods utilized in this work were presented in Table 1. The proposed method does require additional computational time compared to standard methods like Lanczos Tikhonov or other two-step approaches like BPD. Compared to TV-based regularization, the computational time was reduced by $\sim$ 2.5 times, given the performance of the proposed method the additional computational burden is justifiable. The computational complexity of all the methods is ${\mathcal {O}}(n^2)$ not including the SVD complexity of ${\mathcal {O}}(n^4)$ (since SVD of system matrix is computed off-line as one time overhead).

Table 1. Recorded computational time (in seconds) for the results presented in Fig. 4.

View Table

The proposed plug and play prior method for improving the photoacoustic imaging is the simplest of priors that could be utilized or in other terms, we are only utilizing the denoising operator with TV regularization (step-3 in Algorithm-2). One could utilize advanced denoising operators to provide improved performance, the main aim of this work is to introduce the plug and play priors and show that these could provide significant improvement in the reconstruction performance especially low SNR in the raw data. Note that the application of these priors does increase the computational complexity, but in here, it has been cleverly applied with SVD making it dimensionality reduced. Recently deep learning methods have been shown to be very effective for image reconstruction in photoacoustic tomography [23,29–31,73]. The proposed plug and play denoisers can easily be combined with the deep learning based methods for improving the reconstruction. The same was attempted in Ref. [74], where the decoupling was achieved using plug & play models and the prior image denoising operator was modeled using residual deep learning networks. Another deep learning based prior multiple self-similarity net (MSSN) was proposed, based on the recurrent neural network (RNN) with self-similarity matching using multi-head attention mechanism and was shown to give excellent performance for reconstructing three-dimensional (3D) magnetic resonance (MR) images [75]. By applying proposed method slice by slice, the current method can easily be utilized for 3D photoacoustic tomographic image reconstruction, in similar lines to Ref. [31]. Since development of plug & play models based on the proximal gradient method (PGM) utilizes only subset of measurements at every iteration, makes it scalable to very large datasets such as 3D image reconstruction [76].

Earlier there were attempts to utilize additional step of post-processing including deblurring (like BPD and NLM) to provide the robustness to noise in the photoacoustic data [32,77]. The method proposed here utilizes plug and play priors with denoising as a prior operator improving the reconstructed photoacoustic images and can be seen as an iterative method as oppose to a two-step approach. The only significant computational burden for this method is in terms of SVD of system matrix. For a given data-acquisition geometry, the system matrix is typically pre-computed and so is the SVD. The developed algorithms were provided as an open source [78] for enthusiastic users.

6. Conclusion

A singular value based plug and play priors method was proposed for improving photoacoustic imaging and the proposed method was shown to be more robust to noise, making it ideal to be used in noisy data cases. The proposed method was compared with state-of-the-art methods such as Lanczos Tikhonov, basis pursuit deconvolution and total variation based regularization. The improvement obtained in CNR, PC and RMSE is as high as 914.97$\%$, 135.33$\%$ and 67.98$\%$. For in vivo cases the improvement obtained was as high as 8.1 dB with reduction of artifacts. Methods of this type, which provide improved performance irrespective of noise level in the data are universally appealing especially in real-time scenarios. This can further be applied for other image reconstruction problems where the computational complexity is high owing to the large dimensionality of the system matrix.