Photoacoustic tomography (PAT) has the ability to provide optical contrast at ultrasonic resolution deep inside biological tissue.¹ PAT uses a nanosecond laser pulse to irradiate the tissue under investigation; the transient light is absorbed by different tissue chromophores, resulting in a small rise in temperature in the tissue.² The temperature rise generates a pressure wave [also called photoacoustic (PA) waves] due to thermoelastic expansion. The generated PA waves will propagate through the tissue and are measured using broadband ultrasound transducers placed outside the biological tissues.^1,2 These acoustic measurements are used to estimate the initial pressure rise inside the tissue by solving an acoustic inverse problem.² The inversion can be performed using analytical- and model-based methods.^3,4

Model-based reconstruction involves inverting a model matrix that is generated either using impulse response or discretizing the solution of wave equation and is known to perform better than analytical schemes in irregular geometry and limited data scenarios.^3,4 Inversion in limited data scenarios is difficult due to the ill-conditioned nature of the problem.^4,5 Therefore typically, prior statistics about the image is applied in the form of regularization during the inversion.^4,6

Another perspective of regularization lies in its ability to define resolution characteristics in the imaging domain.^7,8 The resolution characteristics are heavily influenced by factors, such as ultrasound transducer sensitivity field, depth-dependent fluence, bandwidth of the detector, and detector position. This work attempts to use the concept of model-based regularization, which is spatially dependent, to mitigate the nonuniform resolution effects arising in PA imaging. The hypothesis is that the nonuniform resolution in the imaging domain can be captured using model resolution. Two regularization schemes within the standard Tikhonov (ST) regularization framework, which use model (system matrix) characteristics, were proposed here. The first one is fidelity-embedded regularization (FER), based on the correlation between the columns of the model matrix, which was proven to be effective in solving inverse problem in electrical impedance tomography.⁹ The second regularization scheme is based on model-resolution matrix, which provides the spatially variant characteristics of the model, that is proven to provide superior results compared to other regularization schemes in diffuse optical tomography.⁸ The performance of the regularization schemes are compared with ST-⁵ and total variation (TV)¹⁰-based schemes using numerical simulations and in-vivo experimental data.

The impulse response of every pixel in the imaging domain was used to build the system matrix to model PA wave propagation as described in Refs. 4 and 11, which was computed by open-source MATLAB k-wave toolbox.¹² The forward model in PAT can be expressed using a linear system of equations⁴ as

The dimension of the system (model) matrix $\mathbf{A}$ is $m\times n$ (here, it is $30720\times 40401$) with $x$ being the column vector containing unknown initial pressure and $b$ being the measured acoustic data.

Estimation of $x$ given a $b$, especially in these limited data cases, requires application of regularization with Tikhonov regularization being the standard one. The cost function (${\mathrm{\Gamma }}_{\mathrm{Tikh}}$) for the Tikhonov regularization scheme is

where $\lambda$ is a regularization parameter providing the balance between residue of the linear equations and expected initial pressure distribution $(x)$. The regularization matrix $\mathbf{R}$ can include prior characteristics of the solution space. A closed-form solution for the minimization of Eq. (2) can be written as Note that the regularization matrix for the zeroth-order Tikhonov (also known as ST) method is an identity matrix, i.e., $\mathbf{R}=\mathbf{I}$.

The proposed FER method incorporates correlations among the column vectors of the system matrix as a regularization matrix during the inversion.⁹ The correlation between the columns of system matrix decreases rapidly as distance between the pixel increases. If two pixels are far from each other, the corresponding columns of system matrix $({\mathbf{A}}_k,{\mathbf{A}}_l)$ are nearly orthogonal. The orthogonal structure gives rise to a weighted backprojected formula that computes the local ensemble average of initial pressure rise at each pixel, which is given as

where $⟨·⟩$ represents the dot product and $x_k^{\mathrm{FER}}$ is the weighted average initial pressure rise at the $k$’th pixel. In a matrix form, to be used in Eq. (3), the $\mathbf{R}$ is a diagonal matrix with $k$’th element being A small modification to the solution by scaling the solution by $\sqrt{1+{\lambda }_{\mathrm{FER}}^2}$ with ${\lambda }_{\mathrm{FER}}$ being the regularization parameter in FER, for ensuring that $x$ being nonzero leads to⁹ Note that, when columns of the system matrix are entirely orthogonal to each other, then FER method will result in ST formulation and when ${\lambda }_{\mathrm{FER}}$ approaches to large value, the solution will be equivalent to filtered backprojection. Thus the ST becomes a special case of FER.

The proposed model-resolution-based regularization (MRR) involves building the model-resolution matrix¹³ from ST case and then deriving a new form of regularization to mitigate the resolution variations arising in the imaging domain. It is defined as^8,13

Each row of $\mathbf{M}$ provides the blur (point spread function) associated with corresponding pixel. The closer the $\mathbf{M}$ is to identity matrix leads to improved resolution in the reconstructed image. The MRR is applied via utilizing a diagonal matrix (${\mathbf{R}}_{\mathrm{MRR}}$) with entries being where ${\mathbf{M}}_{ii}$ is the diagonal entries of $M$ [Eq. (7)]. The normalization used in Eq. (8) minimizes the effect of $\lambda$ in Eq. (7) while computing the diagonal values of model-resolution matrix [${({\mathbf{R}}_{\mathrm{MRR}})}_{ii}$] in Eq. (8). Thus the solution in this case will be In an ideal case, ${\mathbf{R}}_{\mathrm{MRR}}$ will be an identity matrix, resulting in ST formulation [given in Eq. (3)]. However, in reality, ${\mathbf{R}}_{\mathrm{MRR}}$ will capture the nonuniformity among spatial locations to provide model-aware regularization. Note that regularization matrices utilized in both FER and MRR are only model (system matrix and/or regularization parameter) dependent and do not depend on projection data or noise level in the data, thus providing regularization purely that is model dependent.

FER and MRR methods are compared with ST and TV regularization methods, the TV method used here was described in Refs. 10 and 14. The efficacy of the different methods described above was quantified by computing root-mean-square error (RMSE) and peak signal-to-noise ratio (PSNR) of the reconstructed images in numerical experiments and signal-to-noise ratio (SNR) in the in-vivo case.¹⁵

Numerical Derenzo phantom (unipolar in nature, having “0” for the background and “1” for the objects) with initial pressure rise of 1 kPa and a realistic numerical breast phantom created from contrast-enhanced magnetic resonance imaging¹⁶ were considered for comparing the quantitative accuracy of reconstruction methods. In-vivo experimental data from mouse brain region are considered to validate the proposed methods. The details about the experimental setup and data acquisition are given in Ref. 17.

Numerical data were generated using a fine grid having a size of $50.1\times 50.1\mathrm{mm}$ (discretized to $1002\times 1002$) pixels, and the imaging region was $20.1\times 20.1\mathrm{mm}$ with a grid size of $402\times 402\text{pixels}$. The acoustic data were generated on a $402\times 402$ grid using an open-source k-wave toolbox.¹² The collected data were reconstructed on a $201\times 201$ grid to avoid inverse crime. Sixty point detectors having 70% bandwidth and a center frequency of 2.25 MHz are placed equidistantly on a circle of radius 22 mm from the center of the imaging region. The data were sampled at 50 ns with the total time samples being 512. The medium was assumed to be acoustically homogeneous, having speed of sound as $1500\mathrm{m}/\mathrm{s}$, with no absorption and dispersion. Gaussian noise was added to the in-silico forward data to result in data SNR being, 20, 30, and 40 dB.

The reconstruction results using different algorithms pertaining to Derenzo phantom [Fig. 1(m)] were presented in Fig. 1. Figure 1(a) shows the reconstruction results using the ST scheme with data having SNR of 40 dB. Figures 1(b) and 1(c) correspond to the reconstruction results with the ST regularization by having an SNR of 30 and 20 dB, respectively. Figures 1(d)–1(f) represent the reconstruction performance of Derenzo phantom with TV-based reconstruction for simulated data having SNR of 40, 30, and 20 dB, respectively. Reconstruction results pertaining to the proposed FER-based formulation for 40-, 30-, and 20-dB SNR in data were provided in Figs. 1(g)–1(i), respectively. Lastly, the reconstruction results with proposed MRR-based reconstructions for numerical data having SNR of 40, 30, and 20 dB were shown in Figs. 1(j)–1(l), respectively. As shown, the performance of MRR [Fig. 1(l)] and FER [Fig. 1(i)] methods in the 20-dB SNR case is superior when compared with the ST- [Fig. 1(c)] and TV-based schemes [Figs. 1(f) and 1(e)]. Spurious background artifacts tend to arise with TV-based reconstruction as indicated by red arrows in Figs. 1(f) and 1(e). Importantly, the proposed MRR scheme is able to reconstruct images with uniform resolution across the entire imaging domain [indicated by red arrow in Fig. 1(j)]; the background has been recovered more accurately with MRR compared to other reconstruction methods. The measured RMSE and PSNR from the different reconstructions using the Derenzo phantom were given in Figs. 1(n) and 1(o), respectively. RMSE and PSNR metrics indicate that the reconstructed images using FER and MRR provide more accurate solutions compared to ST- and TV-based regularization.

Since the Derenzo phantom is unipolar in nature, the efficacy of the different methods using a realistic breast phantom case [Fig. 2(m)] was investigated and the results were presented in Fig. 2 (which are similar to results presented in Fig. 1). The results demonstrate that fibrous regions were reconstructed better using the proposed MRR scheme than the ST- and TV-type regularizations [Figs. 2(c) and 2(f)], where the noise seems to be amplified. Reconstructions pertaining to FER method with simulated data having SNR of 30 and 20 dB [Figs. 2(h) and 2(i)] were seen to have reduced background artifacts compared to other reconstruction methods as indicated by the red arrow showing its ability in providing robust reconstructions. Uniform resolution across the imaging domain was achieved using the MRR scheme [same is indicated by arrows in Figs. 2(j) and 2(k)]. Similar to previous case, it is clear that the proposed spatial varying regularization schemes, i.e., FER and MRR, outperform ST- and TV-based methods. The reconstructed initial pressure values did not reach to the expected value due to various approximations imposed while performing image reconstruction. The speed of sound was kept constant and the medium was assumed to be homogeneous; the detectors used here were point detectors. Along with these approximations, only a limited number of detectors (60 in here) were considered. Note that utilization of large number of tomographic detector positions will result in more uniform resolution characteristics in the image domain and even the reconstruction accuracy will be high. Since the objective of this work was to demonstrate the potential of the proposed MRR algorithm in mitigating the nonuniform resolution characteristics issue arising in limited detector position case, the case of many detector positions was not considered here.

The central row of model-resolution matrix [computed using corresponding $\mathbf{R}$ instead of $\mathbf{I}$ in Eq. (7)] pertaining to ST, FER, and MRR methods was shown in Figs. 3(a)–3(c), respectively. The central pixel (having highest value) in Fig. 3 represents the row being plotted, and rest pixels correspond to the off-diagonal elements. Ideally, one would expect to see an impulse (center pixel value being 1 and rest pixels being 0). From Fig. 3, it is evident that the off-diagonal elements are closer to zero in case of MRR than Tikhonov and FER. Further ${\Vert \mathrm{diag}(\mathbf{I})-\mathrm{diag}({\mathbf{M}}_{\mathrm{res}})\Vert }_2$ was calculated for each method and found to be 199.81, 199.31, and 173.1 for ST, FER, and MRR, respectively (ideal value being zero). The values of the norm indicate that MRR inversion would result in more resolved solution compared to other methods.

Reconstruction using the in-vivo experimental data using ST, TV, FER, and MRR schemes was shown in Figs. 4(a)–4(d), respectively. The structures shown by the red arrow in Figs. 4(c) and 4(d) were well reconstructed. The computed SNR values (in dB) for each method were shown correspondingly below each reconstructed image. From the reconstruction results in Fig. 4 and SNR values, it is evident that using spatially variant regularization, such as FER and MRR, could potentially enable generation of accurate image representation. Note that the SNR was reported in the dB scale, the improvement is atleast a factor of 1.58 in this in-vivo case using the proposed methods. The computational time for FER, MRR, Tikhonov, and TV is 0.41, 0.38, 0.37, and 9.52 s, respectively, by having the system matrix inversion precomputed. The parameters used for the different reconstructions are ${10}^{-2}$ (20 dB), ${10}^{-3}$ (30 dB), and ${10}^{-4}$ (40 dB) for Tikhonov, ${10}^{-3}$ (20 dB), ${10}^{-3}$ (30 dB), and ${10}^{-4}$ (40 dB) for FER, and ${10}^{-2}$ (20 dB), ${10}^{-4}$ (30 dB), ${10}^{-4}$ (40 dB), and $\mu =0.01$ for MRR, which were set to result in best figure of merit values. The reconstruction parameters, i.e., $\lambda$ and $\mu$, will affect the reconstructed image quality. It is also important to note that automatically choosing multiple parameters will tend to be computationally expensive. Therefore, in this work, these parameters were chosen empirically to result in best possible reconstructions. $\lambda$ is the regularization parameter, which provides a balance between the residual and the expected initial pressure distribution. The parameter $\mu$ provides weight to the normalized spatially varying regularization, which is fixed for all noise cases and phantoms. Higher value of $\mu$ will suppress the effect of spatially varying regularization. The value of $\mu$ is kept as low as possible to provide the best reconstruction. Moreover, these parameters largely depend on the absolute values in the system matrix.

Most regularization schemes present in the literature have limitation in terms of handling spatially varying sensitivity of PA imaging arising due to placement of limited acoustic detectors at the boundary of tissue. In this work, for the first time, spatially varying regularization schemes were proposed to address this issue pertaining to nonuniform resolution characteristics arising in the imaging domain. The proposed regularization schemes were deployed in Tikhonov regularization framework with ST being a special case, making them universally appealing. The results presented in this letter indicate that these model-based regularization schemes were able to outperform ST- and TV-based regularization as they were encapsulating the spatially varying sensitivity of the imaging domain. The MRR was found to generate improved reconstructions (see Figs. 1 and 2, PSNR improvement being as high as 44%) and provides homogeneous resolution characteristics across the imaging domain. The performance of proposed FER relied heavily on the detector positions and the independent nature of the data-collection strategy employed, as the detector positioning will reflect the orthogonality in the columns of system matrix. Typically, FER provides accurate solution with less parameter tuning though the PSNR of FER method is less than the MRR method. Therefore, FER is recommended for fast imaging scenarios. When computing time is not a constraint and more accurate reconstruction is desired, then MRR is a preferred reconstruction method. Note that the quantitative results (SNR, PSNR, and RMSE) clearly demonstrate the robust nature of proposed FER compared to other standard reconstruction methods, such as ST or TV, with more than 2 dB improvement in the case of in-vivo mice brain data. The standard model-based reconstruction in PAT is based on TV,¹⁸ this work with spatially varying regularization was found to outperform standard model-based (Tikhonov/TV) reconstruction. The developed code along with data generation was made available as an open source¹⁹ for enthusiastic users.