## Abstract

Photoacoustic tomography is a promising imaging modality offering high ultrasonic resolution with intrinsic optical contrast. However, quantification in photoacoustic imaging is challenging. We present an algorithm for quantitative photoacoustic estimation of optical absorption and diffusion coefficients based on minimizing an error function between measured photoacoustic channel data and a calculated forward model with a multiple-illumination pattern. Unlike many other algorithms, the proposed method does not require the erroneous assumption of ideal tomographic reconstruction of initial pressures and to our knowledge is the first demonstration of the efficacy of multiple-illumination photoacoustic tomography requiring only transducer channel data. Simulations show promise for numerically robust optical property estimation as illustrated by well-conditioned Hessian singular values in 2D examples.

©2012 Optical Society of America

## 1. Introduction

Photoacoustic tomography (PAT) holds significant promise for high-resolution optical tomographic imaging in optically-scattering tissues. Often, two steps are involved in PAT reconstruction problems. First, the absorbed energy (or initial pressure) distribution is estimated from the measured acoustic signals. Second, one seeks to reconstruct optical properties: the optical absorption coefficient and scattering coefficient map, as well as the Grüneisen parameter distribution based on the results of the first step. This is often called quantitative PAT, or quantitation of PAT. The first step is a well-known acoustic inverse problem and thoroughly studied [1–3]. The second step, an optical inverse problem is, however, non-trivial [4]. This is because of several reasons. The local initial pressures generated when absorbed light pulses are converted to acoustic signals are proportional to not only the local optical absorption coefficient, but also the local laser fluence, which is in turn a complex nonlinear function of the distributed optical properties of the medium. This nonlinear inverse problem of estimating optical properties from photoacoustic data is further complicated by potential ill-posedness: a given photoacoustic absorbed energy distribution may be due to non-unique absorption-scattering distribution pairs. Additionally, when the Grüneisen parameter is unknown as a function of space, Bal and Ren argue that based on diffusion theory, with a single wavelength, only two of the three coefficients in quantitative PAT can be reconstructed uniquely [5].

Quite a number of methods were proposed for quantitative PAT. Studies have been conducted in recovering absorption coefficient distributions with scattering background as a priori information [6–9]. Yin et al. [10] iteratively estimated absorbed energy density with PAT and fluence distribution with diffuse optical tomography. Jetzfellner et al. [11] studied the experimental performance of an iterative scheme and concluded that scattering coefficient errors may lead to failure of the method. Cox et al. [12] recovered both absorption and scattering distribution with multiple optical wavelengths when known wavelength dependence of the optical scattering exists. Efforts were also made later on to decouple Grüneisen coefficients in quantitative PAT. Bal and Ren [13] proposed a scheme for reconstructing absorption, diffusion and Grüneisen coefficients distributions simultaneously with measured data by different wavelengths and *a priori* information on the form of the coefficients stably.

Recently we demonstrated a linearized non-iterative algorithm, namely multiple-illumination photoacoustic tomography (MIPAT), for recovering absorption-scattering distributions using multiple illuminations [14,15]. This approach showed that multiple sources can significantly mitigate absorption-scattering non-uniqueness as demonstrated by simulations. The algorithm we reported used a diffusion-regime ratiometric approach and assumed that initial pressures could be reconstructed in an ideal way. Unfortunately, artifacts in the reconstructed initial pressure distribution can lead to undesirable errors when solving for optical properties. Despite these limitations, the technique showed promise for multiple illumination photoacoustic data for more complex problems.

In this paper, we propose an updated version of the MIPAT algorithm that does not require an ideal initial pressure reconstruction, that instead uses channel data from an ultrasound transducer array as the measurements. The proposed algorithm employs an iterative scheme, therefore it is capable of recovering relatively stronger heterogeneities. While we presently restrict ourselves to the diffusion-regime in 2D, the approach can potentially be generalized to include radiative transport models and acoustic propagation in 2 or 3 dimensions.

## 2. Theory

#### 2.1. Light propagation model

We use the diffusion equation of optical transport as the light propagation model, which can be written as [16]

*q*is the photon density source strength and

*c*is light speed in the medium.

*D*is the diffusion coefficient, which is defined as a function of optical properties by $1/3({\mu}_{a}+{\mu}_{s}^{\text{'}})$. Here ${\mu}_{a}$ and ${\mu}_{s}^{\text{'}}$are the optical absorption coefficient and reduced scattering coefficient, respectively. A time-independent form of this equation is often sufficient when laser-pulse widths are significantly longer than the average random-walk time of photons through the tissue, in which case the time-derivative is neglected. To ensure that the diffusion model is valid, we require${\mu}_{s}^{\text{'}}\gg {\mu}_{a}$. With this hypothesis, fluence $\Phi $ will be almost isotropic. To calculate the fluence perturbations due to background heterogeneity at location $r$, we use the relations [16]

#### 2.2. Reconstruction of the optical properties with ultrasound channel data

Consider that we have a photoacoustic tomography system with *M* ultrasound transducers and *S* illumination patterns. If each detector acquires *T* time points at a given sampling frequency, a column vector of observed pressure measurements may be constructed as

*i*, due to source

*k*. The composite index $\left\{\tau ik\right\}$ can also be assigned the index $\rho $ to map uniquely to the ${\rho}^{th}$ element of this vector such that$\rho =\tau +(i-1)T+(k-1)MT$.

We want to compare these observed photoacoustic channel measurements with forward-model computed values ${p}^{c}$ which are a function of the optical parameters $u={[\begin{array}{cc}{\mu}_{a}& D\end{array}]}^{T}$, where ${\mu}_{a}$ and **D** are rastorized column vectors of the optical absorption and diffusion coefficient maps. The objective is to find best estimate$\widehat{u}$ for $u$ such that

*i*due to source illumination

*k*follow from the superposition of propagated initial pressure signals:

*k*. To model the electromechanical response of the transducer, ${f}^{k}({r}_{i},u,{t}_{\tau})$ is convolved with impulse response ${h}_{r,x}\left(t\right)$. Thus in the expression for ${p}_{\left\{\tau ik\right\}}^{c}$ we can use a filtered version of the initial pressure distribution written as ${p}_{o,k}(r,t)=\Gamma {\mu}_{a}(r){\Phi}_{k}(r,u){h}_{r,x}\left(t\right).$ In future work, the computed channel data may also be modified to include finite aperture effects and acoustic attenuation. By considering that ${f}^{k}({r}_{i},u+\delta u,{t}_{\tau})\cong {f}^{k}({r}_{i},u,{t}_{\tau})+{J}_{i}^{k}\delta u$, where $J$ is the Jacobian matrix, and requiring that $\partial \epsilon (u)/\partial \delta u=0$we arrive at the iterative algorithm $\left({J}^{T}J\right)\delta u={J}^{T}\left({p}^{o}-{p}^{c}(u)\right),$ which can be written aswith$\Delta \equiv {p}^{o}-{p}^{c}(u)$. For each iteration of the algorithm, $\delta u$is computed as an update vector. Once the optical properties are updated, new ${p}^{c}$ vectors and$J$ are computed, then a new update vector $\delta u$ is computed, and the iterations continue until successive iterations of the algorithm yield little change. With Hessian matrix approximated as $H={J}^{T}J$, a Gauss-Newton step may be added to improve convergence speed such that $\delta u=\xi {H}^{-1}[{J}^{T}\Delta ]$, where $\xi $ is a scalar chosen via a line-search algorithm to further minimize the error functional for each step. The Jacobian matrix has elements

*k*. It is straightforward to show that

*i*. $\partial {E}_{\ell}^{k}/\partial {u}_{j}$ represents optical absorbed energy variation with respect to optical property perturbations and it can be written as:

Given that the diffusion coefficient may be written as $D=1/3({\mu}_{a}+{\mu}_{s}^{\text{'}})$ and ${\mu}_{s}^{\text{'}}\gg {\mu}_{a}$, we have $D\approx 1/3{\mu}_{s}^{\text{'}}$. Therefore,

*TMS*× 1. For an

*N*×

*N*2D grid of optical properties, $u$and $\delta u$ have dimensions 2N

^{2}× 1. Jacobian matrices have dimensions

*TMS ×*2

*N*

^{2}. These matrices can be quite large. The Hessian

**H**of size 2

*N*

^{2}× 2

*N*

^{2}has quadrant symmetry, reducing memory requirements. With the above notation, the $\rho th$ row of $J$ is a 1 × 2

*N*

^{2}vector given as${J}_{P}={J}_{\left\{\tau ik\right\}}={[{B}^{k}{\alpha}^{\left\{\tau ik\right\}}]}^{T}$where for a given

*k*, ${B}^{k}$is a 2

*N*

^{2}× 2

*N*

^{2}matrix and for a given $\left\{\tau ik\right\}$, ${\alpha}^{\left\{\tau ik\right\}}$ is an 2N

^{2}× 1 column vector. The Jacobian needs not to be stored, but calculated row by row to find $b={J}^{T}\Delta ={\displaystyle {\sum}_{\rho =1}^{TMS}{J}_{\rho}^{T}}{J}_{\rho}$ and $H={\displaystyle {\sum}_{\rho =1}^{TMS}{J}_{\rho}^{T}}{J}_{\rho}$ in a loop over $\rho $ (a loop over the indices$\left\{\tau ik\right\}$). Future Gradient-based optimization methods may circumvent the necessity of using such large matrices.

#### 2.3. Inversion

Rewrite Eq. (7) as

where$b={J}^{T}\Delta $. Now estimating the unknown parameters $\delta u$ is an inverse problem. Different techniques can be employed to solve the problem. For the*i*iteration, we minimize the following cost function

^{th}**W**is a 2-dimensional regularization operator matrix or weighting matrix. Here we use the second order derivatives as a smoothing weighting factor. For each iteration, the least-squares minimum solution is found with$\nabla {J}_{i}=0$, which yields$\delta {\widehat{u}}_{i}={\left({H}_{i-1}^{T}{H}_{i-1}+\lambda {W}^{T}W\right)}^{-1}{H}_{i}^{T}b.$ Solution $\delta {\widehat{u}}_{i}$is then used to update optical fluence estimation for next iteration. Initial guess of the parameters to be estimated is chosen to be the background absorption and scattering distribution.

## 3. Numerical simulation

With the framework above, we simulated multiple-illuminations of 2D absorption and scattering distributions 2 cm × 2 cm in a grid of 30 × 30 pixels of size $\Delta x=\Delta y=667\mu m$. Configuration of the simulation numerical model and the true scattering-absorption distribution is shown in Fig. 1
. The reduced scattering coefficient of the turbid background is taken as 10 cm^{−1} everywhere except a scattering perturbation of 1 cm^{−1}. The background absorption coefficient is taken as ${\mu}_{a}$ = 0.1 cm^{−1}, with two regions of 10% perturbation. We use 4 illumination sources and an array of 64 ultrasound point detectors distributed 1 cm away from the object (16 on each side of the object). The optical sources are placed 3 mm away from the object to validate the light transport model in the diffusion regime. Modeling incident light as isotropic point sources is equivalent to the pencil beam interrogation that may be used experimentally, by the similarity principle [16]. Keeping these sources far from the simulation area ensures that the diffusion approximation will hold. We sample pressure signals generated by photoacoustic effects with a temporal sampling frequency of 15MHz. The transducer electromechanical response ${h}_{rx}(t)$ is modeled as a low-pass Hanning filter of width 7 time samples. 1% normally distributed white noise (zero mean and 1% of the mean of data) is added to measurements. Figure 2
shows some ‘measured data’ we generated with the simulation model. (a) is the optical normalized fluence generated with a point source on top of the grid system. (b) and (c) are the fluence perturbations due to only absorption and scattering with one point optical source on top of the grid system, respectively. (d) is the total fluence perturbations with illumination from the same optical source. (e) illustrates the normalized ultrasound data collected with the 64 detectors employed, but displayed in different forms.

The time vector has maximum length 283 samples, so with 4 sources the Jacobian is of size 72448 × 1800, the pressure vectors are of length 1800, and the Hessian matrix is 1800 × 1800.

Reconstructed optical absorption and scattering distributions are shown in Fig. 3
. Both *D* and *μ _{a}* distributions are faithfully recovered with the proposed method. Generally

*μ*is better recovered than scattering features. This supports conclusions by other researchers [17]. We believe this is due to stronger dependence of photoacoustic signals on optical absorption. This is also confirmed with faster convergence of optical absorption, as is shown in Fig. 4 . A single iteration provides a reasonable first estimate of both

_{a}*μ*and

_{a}*D*. Additional iterations improve absorption distribution estimates in particular, however, diffusion coefficient begins diverging and this becomes appreciable after 10 iterations. The inclusion in the images of the diffusion coefficient initially have an erroneously low value then iterations improve this mean value, however, successive iterations are sensitive to noise amplification resulting in an overall diverging trend. By simulation and experimental work, Jetzfellner et al. [11] found that a fixed point iteration scheme for recovering absorption heterogeneities could lead to divergence with over-iteration. Their work involved a single illumination pattern. They suggested that

*a priori*information on the imaged object may be necessary. Our simulation showcased improved convergence in absorption distributions, albeit with simulated data. Future work is needed to stabilize diffusion coefficient reconstructions to mitigate iteration-to-iteration noise amplification. As a comparison, we also show results with the ratiometric method described previously in [15] in the bottom row. Only data with signal-to-noise ratio above 60 dB were used in this case. As is shown, the error in the reconstructed absorption map is unacceptably high for the ratiometric approach. Reconstruction of absorption failed. This is because the diffusion coefficient D has higher values than that in the previous paper and cross-talk between absorption and scattering begin to appear. Where the proposed method can alleaviate the cross-talk, the ratiometric method cannot. Scattering features are recovered better than absorption. Readers can refer to [15] for discussion.

The singular values of the Hessian are analyzed as an indicator of inversion stability. The first-iteration Hessian with 1, 4 and 8 sources are shown in Fig. 5 . We found that conditioning improves with multiple illumination sources compared to one source alone (condition numbers are shown in Table 1 ). Detector count and positioning is also advantageous for inversion stability. 64 detectors yield better conditioning than circumferentially positioned 16 detectors. We also found that equally distributed sources perform better than densely positioned sources on one side of the object (data not shown). Further iteration does not appear to significantly impact conditioning. Multi-source configuration is commonly used on diffuse optical tomography, therefore we also compared the conditioning of our proposed methods with continuous wave diffuse optical tomography (CW-DOT). With the same number of optical sources and detectors, the proposed method is orders of magnitude less ill-conditioned.

## 4. Conclusions and discussion

Our algorithm appears to be a promising way of reconstructing optical properties quantitatively using photoacoustic channel data. The use of multiple optical sources helps mitigate absorption-scattering non-uniqueness. The proposed algorithm performs image reconstruction and quantitative optical property estimation together in an iterative scheme.

Our approach is presently based on diffusion theory calculations. However, it may be generalized to incorporate the radiative transport equation (RTE). Though RTE outperforms light diffusion models in reconstruction of optical properties in photoacoustic tomography [17], the diffusion approximation is a reasonable starting point. In our simulation setup, we also took into consideration the feasibility of the diffusion approximation and located light sources far from the object to be imaged. The Born approximation is utilized in our work, limiting our method to be validated with relatively small perturbations.

Additional nonlinear optimization methods should be explored to further improve reconstruction quality. Various constraints should be incorporated during inversion. One example is that Gao et al. [18] introduced the Bregman method combined with total variation regularization for reconstruction of optical properties in photoacoustic imaging. The${\ell}_{1}$ total variation regularization preserved sharp edges and therefore piece-wise continuous features in the recovered map. Simulation work demonstrated that their method outperformed the Jacobian matrix-based methods in terms of computational efficacy.

Our work has limitations. We simplified our work by assuming the Gruneisen parameter is constant throughout the object to be imaged. In practical scenarios, this may lead to inaccuracy in reconstruction of optical properties. We used the same theoretical model for both forward and inverse problems, but to avoid the ‘inverse crime’ we did add noise to our ‘measured’ data. It should also be straightforward to incorporate more realistic transducer aperture, directivity, and electromechanical response effects into future experimental work.

## Acknowledgments

We gratefully acknowledge funding from NSERC (355544-2008, 375340-2009, STPGP 396444), Terry-Fox Foundation and the Canadian Cancer Society (TFF 019237, TFF 019240, CCS 2011-700718), the Alberta Cancer Research Institute (ACB 23728), the Canada Foundation for Innovation, Leaders Opportunity Fund (18472), Alberta Advanced Education & Technology, Small Equipment Grants Program (URSI09007SEG), University of Alberta Startup Funds, and China Scholarship Council scholarships for graduate student Peng Shao.

## References and links

**1. **L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron. **14**(1), 171–179 (2008). [CrossRef]

**2. **M. Xu and L. V. Wang, “Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **67**(5 Pt 2), 056605 (2003). [CrossRef] [PubMed]

**3. **Z. Guo, C. Li, L. Song, and L. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt. **15**(2), 021311 (2010). [CrossRef] [PubMed]

**4. **B. Cox, J. G. Laufer, S. R. Arridge, and P. C. Beard, “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt. **17**(6), 061202 (2012). [CrossRef] [PubMed]

**5. **G. Bal and K. Ren, “Multi-source quantitative photoacoustic tomography in a diffusive regime,” Inverse Probl. **27**(7), 075003 (2011). [CrossRef]

**6. **B. T. Cox, S. R. Arridge, K. P. Köstli, and P. C. Beard, “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt. **45**(8), 1866–1875 (2006). [CrossRef] [PubMed]

**7. **Z. Yuan and H. B. Jiang, “Quantitative photoacoustic tomography: recovery of optical absorption coefficient maps of heterogeneous media,” Appl. Phys. (Berl.) **88**, 231101 (2006).

**8. **J. Ripoll and V. Ntziachristos, “Quantitative point source photoacoustic inversion formulas for scattering and absorbing media,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **71**(3 Pt 1), 031912 (2005). [CrossRef] [PubMed]

**9. **B. Banerjee, S. Bagchi, R. M. Vasu, and D. Roy, “Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map,” J. Opt. Soc. Am. A **25**(9), 2347–2356 (2008). [CrossRef] [PubMed]

**10. **L. Yin, Q. Wang, Q. Zhang, and H. B. Jiang, “Tomographic imaging of absolute optical absorption coefficient in turbid media using combined photoacoustic and diffusing light measurements,” Opt. Lett. **32**(17), 2556–2558 (2007). [CrossRef] [PubMed]

**11. **T. Jetzfellner, D. Razansky, A. Rosenthal, R. Schulz, K.-H. Englmeier, and V. Ntziachristos, “Performance of iterative optoacoustic tomography with experimental data,” Appl. Phys. Lett. **95**(1), 013703 (2009). [CrossRef]

**12. **B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A **26**(2), 443–455 (2009). [CrossRef] [PubMed]

**13. **G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inverse Probl. **28**(2), 025010 (2012). [CrossRef]

**14. **R. J. Zemp, “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt. **49**(18), 3566–3572 (2010). [CrossRef] [PubMed]

**15. **P. Shao, B. Cox, and R. J. Zemp, “Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography,” Appl. Opt. **50**(19), 3145–3154 (2011). [CrossRef] [PubMed]

**16. **L. V. Wang and H.-I. Wu, *Biomedical Optics: Principles and Imaging*, (Wiley, 2007).

**17. **B. Cox, T. Tarvainen, and S. Arridge, in *Contemporary Mathematics* (Book News Inc., 2011), pp. 1–12.

**18. **H. Gao, H. Zhang, and S. Osher, “Bregman methods in quantitative photoacoustic tomography,” CAM Report 10–42, (July 2010).