Open Access
Add to BibSonomyAbstract
The model-based image reconstruction approaches in photoacoustic tomography have a distinct advantage compared to traditional analytical methods for cases where limited data is available. These methods typically deploy Tikhonov based regularization scheme to reconstruct the initial pressure from the boundary acoustic data. The model-resolution for these cases represents the blur induced by the regularization scheme. A method that utilizes this blurring model and performs the basis pursuit deconvolution to improve the quantitative accuracy of the reconstructed photoacoustic image is proposed and shown to be superior compared to other traditional methods via three numerical experiments. Moreover, this deconvolution including the building of an approximate blur matrix is achieved via the Lanczos bidagonalization (least-squares QR) making this approach attractive in real-time.
© 2014 Optical Society of America
1. Introduction
Photoacoustic (PA) imaging is an emerging, noninvasive, in vivo biomedical optical imaging modality that combines both optics and ultrasonic physics [1–3]. Typically a nanosecond laser pulse irradiates biological tissue at near-infrared (NIR) wavelength for deep penetration. The rise in temperature (in the order of milli Kelvin) due to the light absorption by tissue causes emission of pressure waves via thermoelastic expansion. The initial pressure rise is proportional to the absorbed optical energy and the Grueneisen parameter (a dimensionless parameter of the tissue). This pressure wave travels in the soft biological tissues as an acoustic wave, also known as PA wave. A wideband ultrasonic transducer (UST) acquires the propagated PA waves outside the tissue boundary. The UST can be either an array of ultrasonic detectors [4, 5] or a single-element ultrasonic transducer [6, 7]. The PA wave that is collected outside the tissue boundary is used to map the initial pressure rise (or the absorbed optical energy density) within the tissue with the help of a reconstruction algorithm. Structural and functional PA imaging for both pre-clinical and clinical applications has been demonstrated in the literature [4, 5, 8, 9]. In addition, with the help of targeted contrast agents, the photoacoustic imaging has been shown to be a strong contender for the molecular imaging [10, 11].
The critical step in PA imaging is the image reconstruction, where the recent emphasis has been on obtaining quantitatively accurate PA images [12–18]. Several photoacoustic image reconstruction algorithms exist in the literature, including analytic algorithms like filtered back projection (FBP), time-reversal, and Fourier transform based reconstruction [12–17]. These algorithms are based on the spherical Radon transform model, which has an inherent limitation of requiring large number of data points around the target object for accurately estimating the initial pressure distribution. Large number of data points in turn require expensive transducer arrays or long data acquisition time (if single element transducers is used). Moreover, analytic reconstruction does not provide the desired quantitative accuracy with less number of data points. To overcome these limitations, several iterative image reconstruction algorithms [15–18] have been proposed to improve the reconstructed image quality, whose requirement for number of measurements are comparatively less.
The typical PA image reconstruction problem can be considered as a source reconstruction problem, similar to the linear inverse problem that is encountered in Positron Emission Tomography (PET). Previous works have used least-square QR (LSQR) based decomposition methods in obtaining accurate reconstruction with less number of detectors [19–22]. Moreover, optimal selection of regularization parameter using LSQR based method provides quantitatively more accurate reconstruction compared to L-curve and GCV based methods [22, 23]. Hence in this work, LSQR based method was deployed to estimate optimal regularization parameter [22].
As with any model-based image reconstruction methods, deployment of regularization parameter will inherently blur the reconstructed images [24, 25]. These reconstructed images can be deblurred for better recovery of the internal tissue structure as long as the source of the blurring is known or can be modeled. Recently, blind deconvolution method has shown to improve the resolution in optical-resolution photoacoustic microscopy (OR-PAM) [26], in which the point spread function (PSF) was estimated. If the PSF is known then a constrained deconvolution method or least square filtering methods could be used for performing the image deconvolution [27]. Compressive sensing techniques like basis pursuit deconvolution (BPD) has been used in the field of image analysis, and found that this mechanism has the ability to retain the structural information [28, 29]. The BPD algorithms deploy l1-norm based minimization of the generalized least-square objective function. Basis pursuit denoising has been used previously in photoacoustic imaging to perform efficient reconstruction using fewer measurements [30, 31].
The main aim of this work is to remove the blur introduced by regularization parameter, where the blurring model can be built by using the model-resolution characteristics, making the image reconstruction problem a two-step approach, with initial step being model-based traditional reconstruction step and second one being the deblurring step. The deblur-ring/deconvolution of the reconstructed image is achieved via BPD. More importantly, the computation of blur matrix is achieved through Least squares QR (LSQR) approach, which provides an computationally efficient dimensionality reduction procedure. Split augmented Lagrangian shrinkage algorithm (SALSA) is a well known algorithm for performing BPD and the same has been used here [32]. The original SALSA algorithm was based on sparsity signal processing toolbox [33]. In this work, it has been shown that model-resolution based deblurring always results in qualitatively/quantitatively more accurate reconstruction compared to traditional model-based reconstruction, provided the regularization parameter was choosen within a reasonable bounds. The proposed method has the distinct advantage of removing the unwarranted bias in the reconstructed PA image introduced by the heuristic choice of regularization parameter.
2. Analytical and LSQR based reconstruction for photoacoustic tomography
2.1. k-wave based time reversal (analytical)
The k-Wave-based time-reversal is a one-step image reconstruction method, which is computationally efficient, for a planar measurement surface. The mathematical model governing photoacoustic tomography is [34]
Here g(y, t) is the time series pressure data measured at transducer’s location y ∈ B at time t. The pt and are the first and second time derivatives of pressure (p) and c(r) represents the speed of the ultrasound propagation (kept constant in here). The Δr is the Laplacian function with respect to space r. The aim here is to estimate the initial value of f (r) = p(r, t) at t=0, given the measured ultrasound data g(y, t) and c(r).There are many approaches to solve this problem, namely Filtered Back-projection, Fourier Reconstruction and Time reversal methods (all of these are analytical in nature) [17]. The k-wave tool box, which is an open-source based on Matlab, provides analytical tools for performing the initial pressure reconstruction [35]. It is important to note that k-wave can perform full-wave and time-reversal based image reconstruction [18, 35], and in this work time-reversal based reconstruction was used for comparison. The time-reversal method states that for a time point L (the longest time for the PA wave to traverse the domain ω), the solution p(r, t) vanishes inside ω for any t > L [34]. A zero initial condition can now be imposed at t = L and boundary conditions equal to the measured data g and solve the above mathematical model in the time reversal direction, thus arriving at t = 0 to the function f (r) = p(r, 0) [34]. In this work, a interpolated measurement vector is used to estimate the initial pressure of the imaging domain under consideration [35].
2.2. System matrix approach
The process of collection of PA waves outside the boundary can be represented as a time-varying causal system [17]. An image having a dimension of n × n is vectorized by stacking the columns into n2 × 1 vector, represented as x. Each column of the system matrix (A, having dimension of m × n2) is the system’s response to a corresponding entry in the vectorized image (x). It is important to note that the columns of the data (time-varying) are also stacked to result in a long vector having a dimension of m × 1. The system matrix is built as a block circulant matrix, as explained in Ref. [22], with limited transducer bandwidth.
The system response for a corner pixel is recorded using k-Wave MATLAB toolbox [35], similar to Ref. [22]. This system response was found using the k-space pseudo spectral solution method with two coupled first order equations [35]. The simulation geometry of data-collection is shown in Fig. 1. A computational grid indicated by a black square in Fig. 1, having a size of 501 × 501 pixels (0.1 mm/pixel) was used and the detectors (small circles) were placed on a circle of 22 mm radius. Sixty detectors were placed equidistantly around this full circle. The detectors were considered to be a point detector having a center frequency of 2.25 MHz and 70% bandwidth, which were given as set of input parameters into the k-wave simulation for the detection mechanism. Although in practice, large area detectors are used, for simplicity, we have assumed the detectors to be point detector. However, the center frequency and the bandwidth of the point detectors are kept same as the large transducers used in practical PAT systems. A perfectly matched layer (PML) was used to satisfy the boundary conditions. The imaging region (filled green square) was restricted to 201 × 201 pixels located at the center. For data collection, the time step was chosen to be 50 ns with recorded number of time steps being 500. For simplicity, the simulations assumed the speed of sound to be 1500 m/s and the medium to be homogeneous with no absorption or dispersion of sound.
Fig. 1 Schematic of Photoacoustic data acquisition set-up with shaded square region indicating the imaging domain.
The forward model of PA imaging can be summarized as
where x is a long column vector having a dimension of n2 × 1 (n = 201) and b is a measurement vector with a dimension of m × 1 (m = 60 × 500). Back projection (analytical) type image reconstruction scheme will be [36] where, T indicates the transpose of the matrix. Backprojection methods are non-iterative in nature making it computationally efficient, but known to provide only qualitative results, and hence may not be optimal for quantitative imaging [36].2.3. Model based method - Tikhonov regularization (ℓ2-norm based)
The model-based reconstruction relies on minimizing the residue of the linear equations along with a regularization functional having a smoothness constraint, known as Tikhonov Regularization, where the objective function can be written as
where, λ is the regularization parameter, providing the balance between residue of the linear equations (first term on the right-hand side) and expected initial pressure distribution (x). The ℓ2-norm is represented by . The function Ω is minimized with respect to x, resulting in, The regularization parameter dictates the reconstructed initial pressure distribution characteristics. Higher regularization tends to over-smooth the image while lower λ values amplifies the noise in the images.2.4. LSQR-based reconstruction method
The LSQR type method was previously used in photoacoustic tomography [22] for estimation of optimal regularization parameter. The LSQR type algorithm for estimation of approximate model-resolution matrix (blur matrix) and then performing the deblurring is the main contribution of this work. The LSQR based dimensionality reduction of A is achieved by using Lanczos bidiagonalization of A as given in Ref. [37]. The left and right Lanczos matrices and the bidiagonal matrix is related to system matrix (A) as [22, 23]
here, Bk represents the lower bidiagonal matrix, Uk and Vk are the left and right orthogonal Lanczos matrices, respectively. β0 is defined as ℓ2-norm of b. The dimensions of left and right orthogonal Lanczos matrices are (m × k) and (n2 × k), with k representing the number of iterations the bidiagonalization is performed. The unit vector of dimension k × 1 is represented by ek (=1 at the kth row and 0 elsewhere). The structure of Uk, Vk are given by [23] and Bk is a bidiagonal matrix having α1,...αk in the main diagonal and β1,...βk is the lower sub-diagonal of the matrix having a dimension of ((k + 1) × k).Using Eqs. 6, 7, and 8 in the argument of Eq. 4, results in [23]
with x(k) being the dimensionality reduced version of x with k < n2 [38]. Substitution of Eq. 10 in Eq. 4 and using the property Uk+1T Uk+1 = I results in, Considering the first order condition of the Eq. 11, the new update equation becomes [22, 23] where is the estimated version of x(k) for a given regularization parameter λ and fixed k as explained in Ref. [38]. For the case, k = n2, xLSQR = xTikh for a given value of λ. For all other cases, xLSQR becomes an approximation to xTikh. It is well-known that there exists an optimal k with k ≪ n2, for which solutions xLSQR and xTikh are close within the precision limits of the computation [23, 37].2.5. Estimation of optimal λ using a LSQR-type method
The LSQR type algorithm can be used for calculating the (Eq. 12) and then multiplying it by the right Lanczos matrix (Vk) to obtain the approximate solution. This kind of evaluation of update turns out to be computationally more efficient, compared to traditional way of finding update using Eq. 5, as dimensionality reduction is achieved using Lanczos algorithm. Hence, a method of estimating the optimal regularization parameter (λ) was previously proposed in Ref. [22, 23] by minimizing the residue of the linear equations i.e., ||AxLSQR − b||2. The estimation of the number of Lanczos iterations was also performed automatically in Ref. [22, 23]. In this work, the optimal number of Lanczos iterations was found to be 25. The Lanczos bidiagonalization was performed using the Matlab based regularization tools [39].
3. Model-resolution based deconvolution
3.1. LSQR based model resolution matrix
The model resolution matrix can be estimated using the bidiagonal matrix. Substituting the Eq. 6 in Eq. 12 results in,
Substituting Eq. 2 in the above equation gives, Using Eq. 10 converts Eq. 14 into, Rewriting Eq. 7 as, Using the above (Eq. 16) in Eq. 15 results in, which relates the expected x(k) with its estimated version ( ). If the regularization parameter (λ) is equal to zero, then only . The regularization parameter is always greater than 0, indicating that . Rewriting Eq. 17 as, with representing the model resolution matrix (also known as blur matrix) for x(k), having dimension of k × k. Note that the aim of our deconvolution/deblurring is to get an estimate of x(k) using M and . If the regularization parameter is equal to 0, then estimating the blur matrix (M) will not be possible (BkT Bk is a rank-deficient highly ill-conditioned matrix).3.2. Basis pursuit deconvolution in photoacoustic tomography
Deconvolution is a procedure of obtaining a close estimate of the unblurred image from the blurred version. There are several approaches proposed for achieving this [40], with basis pursuit deconvolution being the state of art [28, 29]. In this approach the penalty function is based on ℓ1-norm, which promotes sparseness and sharp features [28,29]. It is known that the regularized reconstructed images are blurry [24], also shown in the previous section (refer to Eq. 18), and deblurring these reconstructed images can be performed by imposing a sparsity constraint. The objective function in this case becomes
where M is the model resolution matrix derived using the least squares QR (LSQR) algorithm (Eq. 19), and (which minimizes Eq. 20 for given λl1) is the deblurred estimate of the . The final estimate from here can be written as where x̃LSQR is the final output (estimate of deblurred xLSQR). It is possible to get exact model-resolution matrix using the Tikhonov formulation (which is computationally expensive) [25], here the estimated model resolution matrix is using LSQR approach which is computationally efficient. It is important to note that is a closer estimate of expected x(k) compared to . The objective function in Eq. 20 can be minimized using Spilt augmented lagrangian shrinkage algorithm (SALSA) [32, 41]. The SALSA algorithm is known to have high convergence speed, when the system matrix (in here M) dimensionality is small, among all existing ℓ1-norm based algorithms achieved via the usage of variable splitting in the minimization problem. This method uses an alternating direction method of multipliers (ADMM), which is based on the augmented Lagrangian method (ALM) [32, 41]. The important steps of SALSA algorithm are given in Algorithm 1. The regularization parameter in the SALSA algorithm (λl1) is chosen as 0.00001 for all cases. The ADMM parameter (d), which has similar size of , is initialized as zero vector. The reconstruction parameter α in the SALSA algorithm weighs the ℓ2-norm of the unknown parameter, which was chosen to be 0.1 for SNR of the data being 40 dB and 10 for all other SNR values. The number of iterations Nit is set to 5000. In here, Nit was chosen to be a large value beyond which the solution did not improve. Note that the computational time for each of these sub-iterations is in the order of milli-seconds. The M and are found using Eq. 19 and Eq. 12 respectively. The initial estimate of is obtained using backprojection-type operation as shown in step-1 of Algorithm-1 (similar to Eq. 2). Soft threshold operation in step-2 indicates an approximation to the ℓ1-norm of and is always non-zero (this is done to avoid the ℓ1-norm penalty from being null, when is close to zero). The soft threshold represents the maximum absolute value among and (0.5λl1)/α. The step-3 in here is a direct translation of Tikhonov-kind of update, which utilizes the normal equations. The step-4 updates the ADMM parameter, ideally when converged to solution there is no update in the d [32, 41]. A flowchart indicating important steps performed in the proposed scheme is given in Fig. 2 for easy understanding, wherein the deconvolution step is clearly indicated.
Algorithm 1:. Spilt augmented lagrangian shrinkage algorithm (SALSA) [41]
4. Figures of merit
4.1. Pearson Correlation (PC)
The performance of the proposed method in terms of improving the reconstructed image quality was evaluated using a quantitative metric, namely Pearson Correlation (PC). The Pearson Correlation is defined as [42]
where x is the expected initial pressure distribution and xrecon is the reconstructed initial pressure distribution using one of the above explained methods. Here σ indicates the standard deviation and COV is the covariance. This measures the degree of correlation between the target (expected) and the reconstructed image, having values ranging from −1 to 1. The higher the value of PC, the better is the detectability of the targets in comparison to the expected image [42].4.2. Contrast to Noise Ratio (CNR)
The other figure of merit that is used is the contrast to noise ratio (CNR), which is defined as,
where μ and σ terms represents the mean and the variance corresponding to the region of interest (ROI, in here the objects) and the background in the reconstructed initial pressure (x). The and represents the area ratio. Higher the value of the CNR better is the reconstruction performance.5. Numerical experiments
It is important to note that measuring actual initial pressure rise in experimental phantom is challenging, leading to difficulties in comparing the quantitative accuracies of the reconstruction methods, making numerical experiments ideal for such comparisons. In order to show the effectiveness of the proposed method, a variation of the derenzo phantom was chosen. This phantom is shown in Fig. 3(a), consists of small and large targets distributions over the imaging region. The targets are grouped as 1–6 according to the sizes as shown in Fig. 3(a). The maximum initial pressure was kept at 1 kPa for the objects.
Fig. 3 (a) Derenzo phantom that is used for the study of resolution characteristics (dimensions are in millimeters). (b–h) Reconstructed photoacoustic images using k-wave interpolated, LSQR with optimal choice of λ, LSQR with heuristic choice of λ, Basis Pursuit Deconvolution (BPD) with optimal choice of λ in LSQR framework, BPD with heuristic choice of λ and 40 dB noise in LSQR framework, BPD with heuristic choice of λ and 30 dB noise in LSQR framework, BPD with heuristic choice of λ and 20 dB noise in LSQR framework, respectively.
Due to the high intrinsic optical contrast (e.g. blood), PA imaging is widely used for visualizing internal blood vessel structures, both in the brain as well as other areas. Another numerical blood vessel phantom as shown in Fig. 4(a), with initial pressure rise as 1 kPa was also used to demonstrate the performance of the algorithm. In terms of comparing the performance of the proposed method with others in terms of reconstruction of sharp edges, a numerical experiment consisting of targets as alphabets ‘PAT’ (shown in Fig. 5(a)) was performed. The k-wave toolbox was used to generate the PA signal with 60 detectors in both the experiments. The collected data was then added with 1% noise in all cases, resulting in signal-to-noise-ratio (SNR) of 40 dB. To test the effectiveness of the proposed scheme for increasing noise levels, the data collection for the derenzo phantom was tested with 20 and 30 dB SNR levels. This data was used for performing the reconstruction using traditional LSQR based reconstruction (without deconvolution) and BPD based two-step approaches (refer to Fig. 2). The k-wave based time reversal reconstruction is also performed for comparison, here this was based on interpolation of the detectors as explained in Ref. [35]. A Linux workstation with dual six-core Intel Xeon processor having a speed of 2.66 GHz with 64 GB RAM has been used for all computations performed in this work.
Fig. 4 (a) Numerical blood vessel phantom that is used for the study (dimensions are in millimeters). (b–f) Reconstructed photoacoustic images using k-wave interpolated, LSQR with heuristic choice of λ, LSQR with optimal choice of λ, Basis pursuit deconvolution (BPD) with heuristic choice of λ in LSQR framework, BPD with optimal choice of λ in LSQR framework, respectively. (g) One-dimensional cross-sectional plot for the results presented in (a),(c),(d),(e), and (f) along the dotted line shown in (a).
Fig. 5 (a) Numerical PAT phantom that is used for the study (dimensions are in millimeters). (b–f) Reconstructed photoacoustic images using k-wave interpolated, LSQR with heuristic choice of λ, LSQR with optimal choice of λ, Basis pursuit deconvolution (BPD) with heuristic choice of λ in LSQR framework, BPD with optimal choice of λ in LSQR framework, respectively. (g) One-dimensional cross-sectional plot for the results presented in (a),(c),(d),(e), and (f) along the dotted line shown in (a).
6. Results and discussion
The reconstructed initial pressure distribution using k-wave time-reversal reconstruction algorithm for derenzo phantom and blood vessel network is shown in Fig. 3(b) and Fig 4(b), respectively. The reconstruction was also performed using the model based reconstruction technique by using the regularization parameter as 0.01 and optimally calculating the regularization parameter using LSQR approach [22]. The reconstruction using the LSQR methods with heuristic choice of regularization for the derenzo and blood vessel phantom is shown in Fig. 3(d) and 4(c), respectively. Optimal choice of regularization parameter was estimated as explained in Ref. [22] and the pressure distribution for this case with derenzo and blood vessel phantom is indicated in Fig. 3(c) and 4(d), respectively. The optimal λ was found to be 0.004 and 0.0036 for the derenzo and the blood vessel phantom, respectively.
The LSQR based approximate model resolution matrix was computed with heuristic choice and optimal choice of regularization parameter. This model resolution matrix was then used to perform the deblurring operation. The initial pressure distribution after performing the deconvolution operation (heuristic choice of λ = 0.01) with the derenzo and the blood vessel phantom is shown in the Fig. 3(f) and 4(e), respectively. The pressure distribution with optimal choice of regularization and deconvolution procedure for the derenzo and the blood vessel phantom is shown in Fig. 3(e) and 4(f), respectively. The initial pressure distribution using the proposed method for varying data SNR levels i.e. 30 dB and 20 dB using the derenzo phantom is shown in Fig. 3(g) and 3(h) respectively. The quantitative metrics (figures of merit) of the reconstructions obtained using other techniques is compiled in Table-1 for comparison (actual reconstructed images are not shown). The results show that the proposed method provides the required noise-tolerance, making it deployable in an experimental scenario, and is able to obtain reasonably more accurate results compared to other methods, when the noise level is high (SNR = 20 dB).
Table 1. The Pearson Correlation (PC) and Contrast to Noise Ratio (CNR) of the reconstructed initial pressure distribution for the results presented in this work.
From Fig. 3(b) it is clearly evident that for limited data points (60 transducer positions, whereas typical PAT system uses 120–250 or more transducer positions [5–8]) k-wave based time reversal reconstruction is inadequate in recovering the original object, even with interpolated data. Also the resolution of the reconstructed objects are position dependent. For group 1 objects the resolution near the scanning centre is much better than those located near the boundary as shown in Fig. 3(b), this degradation in resolution is due to the low numbers of detector position used for data collection as well as the interpolation used. However, model based reconstruction was able to produce location independent object resolution as seen in Fig. 3(c–h).
Another challenge in photoacoustic image reconstruction is the object shape recovery. For object groups 3–6, which are large in size, the reconstructed object looks like donut/ring [Fig. 3(b–d)]. This can be attributed to the limited bandwidth of the detection system. As described earlier, the transducer used for the study have 2.25 MHz center frequency with 70% bandwidth. However with BPD based reconstruction, the object shape is recovered [Fig. 3(e–h)] very close to the original shape. It can be seen that the proposed BPD based reconstruction was able to provide better quantitative accuracy compared to analytical (k-wave time-reversal) based and model based methods (LSQR with heuristic λ and LSQR with optimal λ). The initial pressure rise values show the quantitative accuracy of the BPD based reconstruction.
The line profile for the target and the reconstructed pressure distribution using blood vessel phantom is shown in Fig. 4(g). The line profile was taken along the dotted line as shown in Fig. 4(a). The line profile for the k-wave based time reversal reconstructed image is not shown as visually it is evident that k-wave based time reversal reconstructed image was very poor in quality without delineating the blood vessel structures. Arrows in Fig. 4(b) indicates the blood vessels which were not reconstructed by k-wave. This is due to the limited data points (60) used in k-wave reconstruction. From the line profile, it is clearly evident that BPD based reconstruction, was able to recover the initial pressure rise upto 0.8 kPa compared to the other methods which gave a quantitative value of only 0.3 kPa. The reconstruction results pertaining to ‘PAT’ phantom are as shown in Fig. 5(b–f). It can be seen that the reconstruction using the proposed method yields better contrast compared to the established methods. The optimal choice of regularization parameter (λ) in this case was found to be 1.205e-4. A line profile is shown in Fig. 5(g) that was taken along the dotted line in Fig. 5(a).
Figure 6 shows a plot of residual error ( ) versus the number of Lanczos iterations for the results pertaining to Fig. 4 in terms of choosing optimal k, in here k = 25. Similar behavior was observed for other cases discussed in this work. Note that the relative difference in the residual error (ωk − ωk−1) beyond k = 25 was less than 10−6 (single precision limits). The optimal k (25 for all cases discussed here) is chosen to be the one, when the change in residual error is in single precision limits. Note that the typical computational time taken for each of the method that was discussed in this work is reported in Table-2. This indicates that performing the additional step of BPD has negligible computational burden (2 seconds).
Fig. 6 Plot of the variation of the Residual error and the number of Lanczos iterations for numerical blood vessel phantom (Fig. 4).
Table 2. Typical computational time (in seconds) for reconstructing the initial pressure distribution for the various methods presented in this work. Note that the time taken for building the system matrix (375 seconds) was excluded from this.
The aim of this work is to introduce a basis pursuit deconvolution approach for improving the reconstructed pressure distribution in PAT, where an approximate blur matrix is built using the efficient LSQR approach. The line profile in Fig. 4 indicates that the BPD method results in similar reconstruction with heuristic and optimal choice of regularization parameter (Figs. 4(e) and (f)). The reconstructed pressure distribution in Figs. 3 and 4 indicate that the inclusion of additional step of deconvolution provides better reconstruction than just performing the image reconstruction using the LSQR type approach. The deconvolution process does not require the optimal choice of regularization parameter (λ) as long as it is chosen within the reasonable bounds. These bounds are the values of λ for which
where is a function of λ. In simple words, unless is close to the expected x(k), this relation will not hold good, which essentially means that the BPD methods works well for all cases where is a meaningful result.There were earlier attempts of performing deconvolution of photoacoustic images [26] using blind deconvolution techniques assuming the blur kernel was not known. This work introduced a framework to build the blur matrix based on model-resolution and applied basis pursuit deconvolution to improve the reconstructed pressure distribution. The building of the model resolution matrix was done using a LSQR approach (using bidiagonal matrix) hence providing the computational efficiency than building the model resolution matrix using the original system matrix as explained in Ref. [25].
It is also important to note that the proposed method was able to provide better quantitation compared all other techniques with limited number of detectors (60). Thus, deblurring holds a promise in providing better reconstructed images with less number of detectors, additionally reducing the acquisition time and system cost. The blur that arises due to the laser excitation pulse width or the limited detector bandwidth can also be corrected in the signal domain using a different approach [43]. There, the PA signal itself has been corrected for blurring using Tikhonov regularization. Such methods can also be used on the recorded PA signals and then can be fed to the reconstruction algorithm, where deblurring due to regularization can be corrected using our proposed technique.
The proposed framework in this work has only addressed the problem of blur induced by the regularization, especially in cases where regularization is necessary for providing meaningful results, example being limited data scenarios. The deblurring process has significantly improved the quantitativeness of the results. To make photoacoustic imaging truly quantitative (reconstruction of optical absorption coefficient is essential), more corrections related to incorporation of attenuation of laser light in tissue, non-uniform illumination of light, depth dependent attenuation corrections, limited bandwidth detector response, and laser pulse width, have to be incorporated into the reconstruction scheme. In this work, the limited bandwidth detector response is incorporated into the system model, steps are being undertaken to build a complete model that can represent the realistic photoacoustic imaging scenario. But, this kind of model building for a stand-alone photoacoustic imaging system may not be feasible in all cases, example being tissue absorption being heterogenous (unknown in in-vivo case), which results in non-uniform illumination to the deeper regions, making the model still incomplete.
7. Conclusions
Photoacoustic imaging has been shown to be an invaluable imaging modality both in pre-clinical and clinical settings. The quantitative accuracy of the photoacoustic images can be considerably improved with the usage of model-based reconstruction techniques, especially in cases of limited data. These model-based reconstruction techniques utilizes regularization, typically ℓ2-norm based, to obtain more meaningful reconstruction results, with a caveat of smoothening (blurring) the solutions. This work utilized the least-squares QR-based framework to build a approximate model-resolution (blur) matrix, which in turn got utilized in basis pursuit deconvolution approach to restore the sharp features of the reconstructed photoacoustic image. More importantly, it was shown that the proposed methodology works well when the regularization parameter is chosen within reasonable bounds (rather than at only optimal value), reducing the computational overhead considerably. It was also shown that using the proposed framework, the quantitative accuracy of the reconstructed photoacoustic image has improved by more than 50%.
Acknowledgments
This work is supported by Department of Biotechnology (DBT) Rapid Grant for Young investigator (RGYI) (No: BT/PR6494/GDB/27/415/2012) and DBT Bioengineering Grant (No: BT/PR7994/MED/32/284/2013). JP acknowledges the support by Microsoft Corporation and Microsoft Research India under the Microsoft Research India PhD Fellowship Award and SPIE Optics and Photonics Education Scholarship.
References and links
1. V. E Gusev and A. A Karabutov, Laser Optoacoustics (AIP, 1993).
2. A. A. Karabutov, E. V. Savateeva, and A. A. Oraevsky, “Optoacoustic tomography: New modality of laser diagnostic systems,” Laser Phys. 13, 711–723 (2003).
3. L. H. V. Wang and S. Hu, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science 335, 1458–1462, (2012). [CrossRef] [PubMed]
4. H. F. Zhang, K. Maslov, G. Stoica, and L. H. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotech. 24, 848–851 (2006). [CrossRef]
5. B. Yin, D. Xing, Y. Wang, Y. Zeng, Y. Tan, and Q. Chen, “Fast photoacoustic imaging system based on 320-element linear transducer array,” Phys. Med. Bio. 49, 1339–1346 (2004). [CrossRef]
6. G. Ku and L. H. V. Wang, “Scanning thermoacoustic tomography in biological tissue,” Med. Phys. 27, 1195–1202 (2000). [CrossRef] [PubMed]
7. M. Pramanik, G. Ku, C. H. Li, and L. H. V. Wang, “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys. 35, 2218–2223 (2008). [CrossRef] [PubMed]
8. M. Pramanik and L. H. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt. 14, 054024 (2009). [CrossRef] [PubMed]
9. K. H. Song and L. H. V. Wang, “Noninvasive photoacoustic imaging of the thoracic cavity and the kidney in small and large animals,” Med. Phys. 35, 4524–4529 (2008). [CrossRef] [PubMed]
10. Y. Wang, X. Y. Xie, X. D. Wang, G. Ku, K. L. Gill, D. P. O’Neal, G. Stoica, and L. H. V. Wang, “Photoacoustic Tomography of a Nanoshell Contrast Agent in the in Vivo Rat Brain,” Nano Lett. 4, 1689–1692 (2004). [CrossRef]
11. A. De la Zerda, C. Zavaleta, S. Keren, S. Vaithilingam, S. Bodapati, Z. Liu, J. Levi, B. R. Smith, T. J. Ma, O. Oralkan, Z. Cheng, X. Y. Chen, H. J. Dai, B. T. Khuri-Yakub, and S. S. Gambhir, “Carbon nanotubes as photoacoustic molecular imaging agents in living mice,” Nat. Nanotech. 3, 557–562 (2008). [CrossRef]
12. P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic and photoacoustic tomography,” European J. App. Math. 19, 191–224 (2008).
13. M. Xu and L. H. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Review E 71, 016706 (2005). [CrossRef]
14. K. P. Kostli and P. C. Beard, “Two-dimensional photoacoustic imaging by use of Fourier transform image reconstruction and a detector with an anisotropic response,” App. Opt. 42, 1899–1908 (2003). [CrossRef]
15. P. Ephrat, L. Keenliside, A. Seabrook, F. S. Prato, and J. J. L. Carson, “Three-dimensional photoacoustic imaging by sparse-array detection and iterative image reconstruction,” J. Biomed. Opt. 13, 054052 (2008). [CrossRef] [PubMed]
16. M. A. Anastasio, J. Zhang, X. Pan, Y. Zou, G. Ku, and L. H. V. Wang, “Half-time image reconstruction in thermoacoustic tomography,” IEEE Trans. Med. Imag. 24, 199–210 (2005). [CrossRef]
17. K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Bio. 57, 5399–5423 (2012). [CrossRef]
18. C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media,” IEEE Trans. Med. Imag. 32, 1097–1110 (2013). [CrossRef]
19. X. L. Dean-Ben, V. Ntziachristos, and D. Razansky, “Acceleration of Optoacoustic Model-Based Reconstruction Using Angular Image Discretization,” IEEE Trans. Med. Imag. 31, 1154–1162 (2012). [CrossRef]
20. X. L. Dean-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate Model-Based Reconstruction Algorithm for Three-Dimensional Optoacoustic Tomography,” IEEE Trans. Med. Imag. 31, 1922–1928 (2012). [CrossRef]
21. A. Buehler, A. Rosenthal, T. Jetzfellner, A. Dima, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic inversions with incomplete projection data,” Med. Phys. 38, 1694–1704 (2011). [CrossRef] [PubMed]
22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).
23. J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40, 033101 (2013). [CrossRef] [PubMed]
24. J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized- likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Process. 5, 1346–1358 (1996). [CrossRef]
25. J. Prakash, H. Dehghani, B. W. Pogue, and P. K. Yalavarthy, “Model-Resolution based Basis Pursuit Deconvolution Improves Diffuse Optical Tomographic Imaging,” IEEE Trans. Med. Imag. 2014 (in press). [CrossRef]
26. J. Chen, R. Lin, H. Wang, J. Meng, H. Zheng, and L. Song, “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” Opt. Express 21, 7316–7327 (2013). [CrossRef] [PubMed]
27. B. R. Hunt, “The application of constrained least square estimation to image restoration by digital computer,” IEEE Trans. Comput. C-22, 805–812 (1973). [CrossRef]
28. E. Candes and M. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25, 21–30 (2008). [CrossRef]
29. J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag. 25, 14–20 (2008). [CrossRef]
30. J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imag. 28, 585–594 (2009). [CrossRef]
31. Z. Guo, C. H. Li, L. Song, and L. H. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt. 15, 021311 (2010). [CrossRef] [PubMed]
32. M. Figueiredo, J. Bioucas-Dias, and M. Afonso, “Fast Frame-Based Image Deconvolution Using Variable Splitting and Constrained Optimization,” IEEE Worskhop on Statistical Signal Processing, Cardiff, Wales, 2009.
33. I. Selesnick, “Introduction to Sparsity in Signal Processing [Connexions Web site],” Available at: http://cnx.org/content/m43545/1.3/. (2012).
34. Y. Hristova, P. Kuchment, and L. V. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inv. Problems 24, 055006 (2008). [CrossRef]
35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]
36. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acous. Soc. Am. 112, 1536–1544 (2002). [CrossRef]
37. C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Software 8, 43–71 (1982). [CrossRef]
38. M. E. Kilmer and D. P. OLeary, “Choosing regularization parameters in iterative methods for ill-posed problems,” SIAM J. Matrix Anal. Appl. 22, 1204–1221 (2001). [CrossRef]
39. P. C. Hansen, “Regularization tools version 4.0 for Matlab 7.3,” Numerical Algorithms 46, 189–194 (2007). [CrossRef]
40. P. C. Hansen, J. G. Nagy, and D. P. O. Leary, Deblurring Images: Matrices, Spectra, and Filtering, 1 (SIAM, Philadelphia, 2006). [CrossRef]
41. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Fast Image Recovery Using Variable Splitting and Constrained Optimization,” IEEE Trans. Image Process. 19, 2345–2356 (2010). [CrossRef] [PubMed]
42. J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse Recovery Methods Hold Promise for Diffuse Optical Tomographic Image Reconstruction,” IEEE J. Sel. Top. Quantum Electron. 20, 6800609 (2014). [CrossRef]
43. N. A. Rejesh, H. Pullagurla, and M. Pramanik, “Deconvolution based deblurring of reconstructed images in photoacoustic/thermoacoustic tomography,” J. Opt. Soc. Am. A 30, 1994–2001 (2013). [CrossRef]
