Translator Disclaimer
25 October 2018 Spatially variant regularization based on model resolution and fidelity embedding characteristics improves photoacoustic tomography
Author Affiliations +
Photoacoustic tomography tends to be an ill-conditioned problem with noisy limited data requiring imposition of regularization constraints, such as standard Tikhonov (ST) or total variation (TV), to reconstruct meaningful initial pressure rise distribution from the tomographic acoustic measurements acquired at the boundary of the tissue. However, these regularization schemes do not account for nonuniform sensitivity arising due to limited detector placement at the boundary of tissue as well as other system parameters. For the first time, two regularization schemes were developed within the Tikhonov framework to address these issues in photoacoustic imaging. The model resolution, based on spatially varying regularization, and fidelity-embedded regularization, based on orthogonality between the columns of system matrix, were introduced. These were systematically evaluated with the help of numerical and in-vivo mice data. It was shown that the performance of the proposed spatially varying regularization schemes were superior (with at least 2 dB or 1.58 times improvement in the signal-to-noise ratio) compared to ST-/TV-based regularization schemes.

Photoacoustic tomography (PAT) has the ability to provide optical contrast at ultrasonic resolution deep inside biological tissue.1 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.2 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.2 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.9 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.8 The performance of the regularization schemes are compared with ST-5 and total variation (TV)10-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.12 The forward model in PAT can be expressed using a linear system of equations4 as

Eq. (1)

The dimension of the system (model) matrix A is m×n (here, it is 30720×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 (ΓTikh) for the Tikhonov regularization scheme is

Eq. (2)

where λ is a regularization parameter providing the balance between residue of the linear equations and expected initial pressure distribution (x). The regularization matrix R can include prior characteristics of the solution space. A closed-form solution for the minimization of Eq. (2) can be written as

Eq. (3)

Note that the regularization matrix for the zeroth-order Tikhonov (also known as ST) method is an identity matrix, i.e., R=I.

The proposed FER method incorporates correlations among the column vectors of the system matrix as a regularization matrix during the inversion.9 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 (Ak,Al) 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

Eq. (4)

where · represents the dot product and xkFER is the weighted average initial pressure rise at the k’th pixel. In a matrix form, to be used in Eq. (3), the R is a diagonal matrix with k’th element being

Eq. (5)

A small modification to the solution by scaling the solution by 1+λFER2 with λFER being the regularization parameter in FER, for ensuring that x being nonzero leads to9

Eq. (6)

Note that, when columns of the system matrix are entirely orthogonal to each other, then FER method will result in ST formulation and when λ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 matrix13 from ST case and then deriving a new form of regularization to mitigate the resolution variations arising in the imaging domain. It is defined as8,13

Eq. (7)

Each row of M provides the blur (point spread function) associated with corresponding pixel. The closer the M is to identity matrix leads to improved resolution in the reconstructed image. The MRR is applied via utilizing a diagonal matrix (RMRR) with entries being

Eq. (8)

(RMRR)ii=Miimax[diag(M)],for  i=1,2,,n,
where Mii is the diagonal entries of M [Eq. (7)]. The normalization used in Eq. (8) minimizes the effect of λ in Eq. (7) while computing the diagonal values of model-resolution matrix [(RMRR)ii] in Eq. (8). Thus the solution in this case will be

Eq. (9)

In an ideal case, RMRR will be an identity matrix, resulting in ST formulation [given in Eq. (3)]. However, in reality, RMRR 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.15

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 imaging16 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×50.1  mm (discretized to 1002×1002) pixels, and the imaging region was 20.1×20.1  mm with a grid size of 402×402  pixels. The acoustic data were generated on a 402×402 grid using an open-source k-wave toolbox.12 The collected data were reconstructed on a 201×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  m/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.

Fig. 1

Comparison of the performance of different reconstruction methods using Derenzo (unipolar) phantom. Images reconstructed with ST regularization scheme are given in the first column with SNR of data being (a) 40, (b) 30, and (c) 20 dB, TV regularization scheme in second column with SNR of data being (d) 40, (e) 30, and (f) 20 dB. Correspondingly, the proposed schemes, FER and MRR, results are presented in third and fourth columns, respectively. Original Derenzo phantom is shown in (m). RMSE and PSNR (in dB) pertaining to results presented from (a) to (l) are shown in (n) and (o), respectively.


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.

Fig. 2

Same effort as Fig. 1 using realistic breast phantom given in (m).


The central row of model-resolution matrix [computed using corresponding R instead of 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 diag(I)diag(Mres)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.

Fig. 3

Row corresponding to central pixel in the model-resolution matrix for (a) ST, (b) FER, and (c) MRR.


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 102 (20 dB), 103 (30 dB), and 104 (40 dB) for Tikhonov, 103 (20 dB), 103 (30 dB), and 104 (40 dB) for FER, and 102 (20 dB), 104 (30 dB), 104 (40 dB), and μ=0.01 for MRR, which were set to result in best figure of merit values. The reconstruction parameters, i.e., λ and μ, 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. λ is the regularization parameter, which provides a balance between the residual and the expected initial pressure distribution. The parameter μ provides weight to the normalized spatially varying regularization, which is fixed for all noise cases and phantoms. Higher value of μ will suppress the effect of spatially varying regularization. The value of μ is kept as low as possible to provide the best reconstruction. Moreover, these parameters largely depend on the absolute values in the system matrix.

Fig. 4

Reconstructed images of in-vivo mice brain data using (a) ST, (b) TV, (c) proposed FER, and (d) proposed MRR. The corresponding SNR (in dB) values of reconstructed images are given below.


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,18 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 source19 for enthusiastic users.


No conflicts of interest, financial or otherwise, are declared by the authors.



J. Xia and L. V. Wang, “Small-animal whole-body photoacoustic tomography: a review,” IEEE Trans. Biomed. Eng., 61 (5), 1380 –1389 (2014). IEBEAX 0018-9294 Google Scholar


Y. Zhou et al., “Tutorial on photoacoustic tomography,” J. Biomed. Opt., 21 061007 (2016). JBOPFO 1083-3668 Google Scholar


A. Rosenthal et al., “Acoustic inversion in optoacoustic tomography: a review,” Curr. Med. Imaging Rev., 9 318 –336 (2013). Google Scholar


J. Prakash et al., “Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography,” Biomed. Opt. Express, 5 1363 –1377 (2014). BOEICL 2156-7085 Google Scholar


M. Bhatt et al., “Exponential filtering of singular values improves photoacoustic image reconstruction,” J. Opt. Soc. Am. A, 33 1785 –1792 (2016). JOAOD6 0740-3232 Google Scholar


Y. Han et al., “Three-dimensional optoacoustic reconstruction using fast sparse representation,” Opt. Lett., 42 (5), 979 –982 (2017). OPLEDP 0146-9592 Google Scholar


B. W. Pogue et al., “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt., 38 (13), 2950 –2961 (1999). APOPAI 0003-6935 Google Scholar


S. H. Katamreddy and P. K. Yalavarthy, “Model-resolution based regularization improves near infrared diffuse optical tomography,” J. Opt. Soc. Am. A, 29 (5), 649 –656 (2012). JOAOD6 0740-3232 Google Scholar


K. Lee et al., “A fidelity-embedded regularization method for robust electrical impedance tomography,” IEEE Trans. Med. Imaging, 37 1970 –1977 (2018). ITMID4 0278-0062 Google Scholar


M. V. Afonso et al., “Fast image recovery using variable splitting and constrained optimization,” IEEE Trans. Imaging Process., 19 (9), 2345 –2356 (2010). Google Scholar


C. B. Shaw et al., “Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt., 18 (8), 080501 (2013). JBOPFO 1083-3668 Google Scholar


B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wavefields,” J. Biomed. Opt., 15 021314 (2010). JBOPFO 1083-3668 Google Scholar


J. Prakash et al., “Model-resolution based basis pursuit deconvolution improves diffuse optical tomographic imaging,” IEEE Trans. Med. Imaging, 33 (4), 891 –901 (2014). ITMID4 0278-0062 Google Scholar


A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision, 20 (1), 89 –97 (2004). JIMVEC 0924-9907 Google Scholar


O. Marques, Practical Image and Video Processing using MATLAB, John Wiley & Sons, New Jersey (2011). Google Scholar


Y. Lou et al., “Generation of anatomically realistic numerical phantoms for photoacoustic and ultrasonic breast imaging,” J. Biomed. Opt., 22 (4), 041015 (2017). JBOPFO 1083-3668 Google Scholar


Y. Jiang et al., “Broadband absorbing semiconducting polymer nanoparticles for photoacoustic imaging in second near-infrared window,” Nano Lett., 17 (8), 4964 –4969 (2017). NALEFD 1530-6984 Google Scholar


C. Zhang, Y. Zhang and Y. Wang, “A photoacoustic image reconstruction method using total variation and nonconvex optimization,” Biomed. Eng. OnLine, 13 117 (2014). Google Scholar


D. Sanny et al., “Spatially varying regularization for photoacoustic tomography,” (2018) September ). 2018). Google Scholar
© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Dween Rabius Sanny, Jaya Prakash, Sandeep Kumar Kalva, Manojit Pramanik, and Phaneendra K. Yalavarthy "Spatially variant regularization based on model resolution and fidelity embedding characteristics improves photoacoustic tomography," Journal of Biomedical Optics 23(10), 100502 (25 October 2018).
Received: 3 August 2018; Accepted: 27 September 2018; Published: 25 October 2018

Back to Top