Photoacoustic (PA) imaging is an emerging, noninvasive, in vivo biomedical imaging modality.^1,2 A nanosecond laser pulse is generally used to irradiate biological tissue, resulting in a temperature rise from optical absorption and subsequently producing pressure waves due to thermoelastic expansion.² The pressure waves are then acquired using a wide-band ultrasonic transducer at various locations around the surface of the tissue. A reconstruction algorithm is deployed that maps the initial pressure rise (proportional to the absorbed optical energy density) within the tissue from the recorded PA signals.³

Several PA image reconstruction algorithms were proposed in the literature,^4,5 including analytical algorithms in the form of filtered back projection (BP) or algorithm based on Fourier transform. Their limitations include the requirement of large amount of data and are limited in terms of quantitative estimation.^6–8 To overcome this limitation, various iterative image reconstruction algorithms have been proposed^6–8 to improve the quantative accuracy of the reconstructed images, at the same time being computationally efficient. Moreover, in case of full-view data sets, the least squares QR (LSQR)-based reconstruction scheme was used, which indirectly provides regularized solution with an added advantage of being highly efficient.⁷ The limited-view data is reconstructed using a standard Tikhonov regularization, which is time consuming and requires an explicit regularization parameter.^5–8

In this letter, we propose a Tikhonov regularization framework based on LSQR decomposition, where $Q$ and $R$ represent an orthogonal and upper triangular matrices, respectively, which uses Lanczos bidiagonalization to provide dimensionality reduction to the system of equations in the case of limited-view data set. This is further used to carry out a simplex method-based optimization procedure to find the optimal regularization parameter. The performance of the proposed method is compared with the generalized cross validation (GCV)⁹ and L-curve⁹ methods, along with the analytical methods, such as BP⁴ and $k$-wave-based time-reversal reconstruction¹⁰ using a numerical blood vessel network phantom.

The system matrix approach has been adopted here to describe the PA data collection process, which can be represented by a Toeplitz matrix of a time-varying causal system. The image (dimension of $n\times n$) is converted into a long vector by stacking all columns one below the other, represented by $x$ (dimension of $n^2\times 1$). The system matrix ($\mathbf{A}$) has a dimension of $m\times n^2$. Here, each column of $\mathbf{A}$ represents the impulse response corresponding to each pixel in the image. Moreover, the time-varying data is stacked to result in a long vector having dimensions $m\times 1$ (which makes the number of rows of $\mathbf{A}$ to be $m$). In order to improve the computation time for building the system matrix, the system response was measured only once for the corner pixel [making $x(\text{corner pixel})=1$ and the rest of the entries were made zero], which forms the first column entry of $\mathbf{A}$. The rest of the columns ($n^2-1$) are filled by using shifting and attenuation properties of the PA signal. This approach assumes that the medium has homogeneous ultrasound properties.

The system response for the corner pixel is recorded using $k$-wave matlab toolbox,¹⁰ which simulates the PA wave prorogation in two dimension. The simulation geometry had a computational grid of $701\times 701\text{pixels}$ ($0.1\mathrm{mm}/\text{pixel}$). Forty detectors were placed in a circular fashion of 34-mm radius. Each detector was assumed to be a point detector with a frequency response of 2.25 MHz as center frequency and 70% bandwidth. The imaging region was restricted to $201\times 201\text{pixels}$ located at the center, resulting in $n^2=40,401$. A time step of 50 ns having 1000 time steps was used in recording each signal (making $m=40,000$). The simulations assumed a sound speed of $1500\mathrm{m}/\mathrm{s}$.

In summary, the forward model of PA imaging can be written as

where $x$ is a long column vector (unknown, representing the initial pressure rise, $p_0$) and $b$ is a measurement vector. The simple BP (analytical) image reconstruction scheme becomes⁴ where $T$ represents the transpose of the matrix. As it is noniterative, this method is computationally efficient but known to provide only qualitative results.⁴

Both BP and time-reversal methods are analytical in nature but lack the quantitative nature of the results.⁵ In cases of limited data, typically a model-based reconstruction is employed which relies on minimizing the data-model misfit along with a regularization function, therefore the objective (cost) function in this case can be written as

where $\lambda$ is the regularization parameter. The ${\ell }_2$-norm is represented by ${\Vert .\Vert }_2^2$. The function $\mathrm{\Omega }$ is minimized with respect to $x$, leading to a direct solution

The GCV method⁹ is the most popular automated approach for estimating the optimal regularization parameter ${\lambda }_{\mathrm{opt}}$ using following function

where the SVD of $\mathbf{A}=\mathbf{H}\mathrm{\Sigma }{\mathbf{G}}^T$ where, $\mathrm{\Sigma }$ is a diagonal matrix containing singular values ($\sigma$). The left and right orthogonal matrices are given by $\mathbf{H}$ and $\mathbf{G}$, respectively. The L-curve is another popular scheme for estimating the optimal regularization parameter.⁹ The corner of L-curve gives the optimal regularization value, and as with GCV method, it does not require any prior information. The LSQR method is one of the variants of the conjugate gradient method used to solve a large system of equations. One of the main contributions of this letter is to use the LSQR-type algorithm to optimally determine the regularization parameter in PA image reconstruction. This is accomplished by using a Lanczos bidiagonalization of the system matrix ($\mathbf{A}$). The left and right Lanczos matrices along with the bidiagonal matrix are related to the system matrix as shown below:^11,12 where $B$ represents the lower bidiagonal matrix, and $U$ and $V$ represent the left and right orthogonal Lanczos matrices, respectively. The unit vector of dimension $k\times 1$ is represented by $e_k$ ($=1$ at the $k$’th row and 0 elsewhere). Note that the dimensions of $U_k$ and $V_k$ are ($m\times k$) and ($n^2\times k$), with $k$ representing the number of iterations in the bidiagonalization procedure. Finally, $u_i$ and $v_i$ represent the left and right Lanczos vectors. The $B_k$ is the bidiagonal matrix having ${\alpha }_1,\dots ,{\alpha }_k$ in the main diagonal and ${\beta }_1,\dots ,{\beta }_k$ in the lower subdiagonal of the matrix with a dimension of [$(k+1)\times k$].

Now the Tikhonov minimization update for the equation for the LSQR-type method¹² is given by [which is equivalent to Eq. (4)]:

Here ${\beta }_0$ is the ${\ell }_2$-norm of $b$. Once, $x^k$ (reduced $x$) is estimated then the initial pressure can be obtained using the relation $p_0=x=V_k{x}^k$.

Determination of the optimal number of Lanczos iterations ($k^{\mathrm{opt}}$) and optimal regularization parameter ${\lambda }_{\mathrm{opt}}$ is given in Algorithm 1. The advantage of LSQR-type method in finding the initial pressure rise distribution $p_0$ lies in its dimensionality reduction capability which makes the update as $x^k$ [Eq. (8)] with $k\ll n^2$. The major role in the entire optimization procedure is characterized by $k$ (number of Lanczos bidiagonalization). This factor determines the size of the bidiagonal matrix, $B_k$ [dimension of $(k+1)\times k$].^11,12 In this work, the optimal number of Lanczos iterations turned out to be 25 obtained using Algorithm 1, making $k^{\mathrm{opt}}=25$. The optimal $\lambda$ is searched within the specified bound (${\lambda }_{\lim }$) and is set to 1000. A gradient-free simplex method type algorithm is used due to its computational efficiency to compute the optimal regularization parameter (${\lambda }_{\mathrm{opt}}$).¹² The ${\lambda }^{\mathrm{opt}}$ is found corresponding to $k=25$ is chosen for LSQR method. Using these in Eq. (8) gave $x^k$ and subsequently $p_0$ ($x$).

Algorithm 1

Algorithm for determining optimal number of Lanczos iterations and optimal regularization.

A Linux workstation with Intel Xeon Dual Quad Core 2.33 GHz processor having 64 GB memory was used in all computations carried out in this work.

In order to show the effectiveness of the proposed method, a numerical blood vessel phantom was chosen. Note that measuring actual $p_0$ in the experimental phantom case is extremely challenging, which makes the comparison of performance of the methods discussed here difficult. Figure 1(a) shows the blood vessel network used as a numerical phantom with a maximum initial pressure rise of 1 kPa. The $k$-wave tool box¹⁰ was used to generate the simulated PA data with 40 detectors around the object of interest. Subsequently, the simulated data had a signal-to-noise-ratio of 40 dB (1% noise) to mimic the real experimental situation. The reconstruction results obtained using various methods are shown in Fig. 1(b)–1(f). The proposed method’s result is shown in Fig. 1(f), with the value obtained for ${\lambda }_{\mathrm{opt}}$ shown in the parenthesis at the top of the image. Figure 1(g) shows the one-dimensional cross-sectional plot for the reconstructed PA image using GCV, L-curve, and LSQR-based methods as well as target image [along the dotted line in Fig. 1(a)], quantitatively showing an improvement of at least 10 times in the recovered $p_0$. As seen from Fig. 1, the performance of LSQR-type method is superior in terms of quantitation compared to their counterparts. The total computational time recorded for all reconstruction methods shown in Fig. 1 are 129.8, 1.3, 7516.5, 7130.7, 444.9 s, respectively (system matrix building time = 181 s). A speed up factor of 17 was achieved by the proposed method compared to GCV method. The results indicate that among the model-based methods, the proposed method (LSQR) is the most efficient and promising technique for real-time imaging.

It is important to note that the LSQR (without explicit regularization) has been extensively used in full-view data cases (where $m\gg n^2$) and shown to be effective in quantitation.^6–8 In limited data cases (where $m\ll n^2$), the Tikhonov minimization scheme is more effective with the important condition of finding an optimal regularization parameter. This work addresses this problem of finding an optimal regularization parameter automatically without prior information with an added advantage of reducing the dimensionality of the problem, making it a highly computationally efficient method.