## Abstract

We describe a novel reconstruction method that allows for quantitative recovery of optical absorption coefficient maps of heterogeneous media using tomographic photoacoustic measurements. Images of optical absorption coefficient are obtained from a diffusion equation based regularized Newton method where the absorbed energy density distribution from conventional photoacoustic tomography serves as the measured field data. We experimentally demonstrate this new method using tissue-mimicking phantom measurements and simulations. The reconstruction results show that the optical absorption coefficient images obtained are quantitative in terms of the shape, size, location and optical property values of the heterogeneities examined.

© 2007 Optical Society of America

## 1. Introduction

Biomedical photoacoustic tomography (PAT) is a potentially powerful imaging method for visualizing the internal structure of soft tissues with excellent spatial resolution and satisfactory imaging depth.^{1–17} While conventional PAT can image tissues with high spatial resolution, it provides only the distribution of absorbed light energy density that is the product of both the *intrinsic* optical absorption coefficient and *extrinsic* optical fluence distribution. Thus the imaging parameter of conventional PAT is clearly not an intrinsic property of tissue. It is well known, however, that it is the tissue absorption coefficient that directly correlates with tissue physiological/functional information. These physiological parameters including hemoglobin concentration, blood oxygenation and water content are critical for accurate diagnostic decision-making.

Several methods reported suggest that it is possible to recover optical property maps when conventional PAT is combined with a light transport model.^{11–14} However, there are several limitations associated with these methods. First, in these methods, one has to know the exact boundary reflection coefficients as well as the exact strength and distribution of incident light source which require careful experimental calibration procedures. It is often difficult to obtain these initial parameters accurately. Second, the recovered results strongly depend on the accuracy of the distribution of absolute absorbed energy density from conventional PAT. As such, to overcome the limitations mentioned above, in this paper we propose a novel reconstruction approach that combines conventional PAT with diffusion equation based regularized Newton method for accurate recovery of optical properties. This work represents the first application of the diffusion equation based iterative nonlinear algorithms that couple the conventional Tikhonov regularization with a priori spatial information-based regularization schemes for reconstruction of absorption coefficient using tomographic photoacoustic measurements. We demonstrate this method using a series of phantom experiments.

## 2. Methods and materials

In our reconstruction method, the absorbed optical energy density is first recovered by a finite element-based PAT reconstruction algorithm.^{10,17} By incorporating the recovered absorbed energy density distribution into the photon diffusion equation, the absorption coefficient map is then extracted using a diffusion equation based regularized Newton method. The core procedure of our PAT algorithm can be described by the following two equations

in which *p* is the pressure wave; *k*
_{0}=*ω/c*
_{0} is the wave number described by the angular frequency, *ω* and the speed of acoustic wave in the medium, *c _{0}; β* is the thermal expansion coefficient;

*C*is the specific heat; Φ is absorbed light energy density that is the product of optical absorption coefficient,

_{p}*µ*and optical fluence or photon density, Ψ (i.e., Φ=

_{a}*µ*Ψ);

_{a}**p**

^{o}=(

*p*

^{0}

_{1},

*p*

^{0}

_{2},…,

*p*

^{0}

*)*

_{M}*,*

^{T}**p**

*=(*

^{c}*p*

_{c}_{1},

*p*

^{c}

_{2},…,

*p*)

^{c}_{M}*, where*

^{T}*p*and

^{o}_{i}*p*are observed and computed complex acoustic field data for i=1, 2…,

^{c}_{i}*M*boundary locations; Δ

*χ*is the update vector for the absorbed optical energy density; ℑ is the Jacobian matrix formed by ∂p/∂Φ at the boundary measurement sites;

*λ*is a Levenberg-Marquardt regularization parameter

^{15}and

**I**is the identity matrix. Thus here the image formation task is to update absorbed energy density distribution via iterative solution of Eqs. (1) and (2) so that a weighted sum of the squared difference between computed and measured acoustic data can be minimized.

To recover the optical absorption coefficient *µ _{a}(r)* from the absorbed energy density, Φ, the photon diffusion equation as well as the Robin boundary conditions (BCs) can be written in consideration of Φ=

*µ*Ψ,

_{a}where *E*(*r*)=1/*µ _{a}*(

*r*),

*D(r)*is the diffusion coefficient which can be written as

*D*=1/(3(

*µ*+

_{a}*µ*′s)) where

*µ*′s is the reduced scattering coefficients,

*α*is the boundary conditions coefficient related to the internal reflection at the boundary, and

*S(r)*is the incident area source term. Note that there is a difference between the BCs used in Eq. (4) and the Robin BCs used in usual diffuse optical tomography: the BCs used here are involved in both the forward and inverse computations, while the Robin BCs employed in diffuse optical tomography are considered only in the forward calculation. Basically there are two methods for solving an elliptic BV problem with Robin BCs. One is based on analytical solution and the other is using numerical methods (e.g., finite element or finite difference). Analytical-based methods are easy to implement.

^{11}However, it is almost impossible to obtain an analytical solution for realistic geometry and optical properties of tissue. Numerical methods are needed for such cases although they are time-consuming. In addition, numerical methods are often combined with iteration-based optimization algorithms.

For the inverse computation, the Tikhonov-regularization sets up a weighted term as well as a penalty term in order to minimize the squared differences between computed and measured absorbed energy density values,^{16}

where **L** is the regularization matrix or filter matrix, *β* the regularization parameter and **E**
_{0} the initial guess of the inverse of optical absorption coefficient. **Φ**
^{o}=(**Φ**
^{o}
_{1},**Φ**
^{o}
_{2},…,**Φ**
^{o}
_{N})^{T} and **Φ**
^{c}=(**Φ**
^{c}
_{1},**Φ**
^{c}
_{2},…,**Φ**
^{c}
_{N})^{T}, where **Φ**
^{0}
_{i} is the normalized absorbed energy density obtained from PAT, and **Φ**
^{c}
_{i} is the absorbed energy density computed from Eqs. (3) and (4) for *i*=1, 2…, N locations within the entire PAT reconstruction domain. It should be noted that the reconstruction of the inverse of optical absorption coefficient using Eqs. (3) and (4) will make the inverse computation easier. The initial estimate of the inverse of absorption coefficient can be updated based on iterative Newton method as follows,

where **J** is the Jacobian matrix formed by ∂Φ/∂E inside the whole reconstruction domain including the boundary zone. The practical update equation resulting from Eq. (6) is utilized with *β*=1,

In addition to the usual Tikhonov regularization, the PAT image (absorbed energy density map) is used both as input data and as prior structural information to regularize the solution so that the ill-posedness associated with such inversion can be reduced. In our reconstruction scheme, we first segment the PAT image into different regions according to different heterogeneities or tissues types using commercial software. We then employ both the distribution of absorbed energy density in the entire imaging domain and segmented prior structural information for optical inversion. The segmented prior spatial information can be incorporated into the iterative process using the regularization filter matrix, **L** shown in Eq. (7). In this study, Laplacian-type filter matrix is employed and constructed according to the region or tissue type it is associated based on derived priors. This filter matrix is able to relax the smoothness constraints at the interface between different regions or tissues, in directions normal to their common boundary so that the co-variance of nodes within a region is basically realized. As such, the elements of matrix **L**, *L _{ij}*, is specified as follows

^{16}:

where *NN* is the total node number within one region or tissue. It should be noted the last term in Eq. (6) is not routinely used in the reconstruction and including the term would reduce the sharpness of known edges given a homogeneous initial guess. Thus the absorption coefficient distribution is reconstructed through the iterative procedures described by Eqs. (3) and (7).

The image formation process described above is tested first using simulated data. The test geometry is shown in Fig. 1a where a two-dimensional (2D) circular background region (50.8 mm in diameter) contained four circular targets (5.08 mm in diameter each). The optical properties for the targets were *µ _{a}*=0.04 mm

^{-1}and

*µ′*=1.0 mm

_{s}^{-1}while optical properties for the background were

*µ*=0.01 mm

_{a}^{-1}and

*µ′*=1.0 mm

_{s}^{-1}. In the simulation, a homogeneous distributed area source is utilized to illuminate the whole imaging domain from its top surface, which is the same as in our experiments (see below). A total of 120 ultrasound receivers were equally distributed along the boundary of background region. While PAT signals carry a wide range of acoustic frequencies, only 50 frequencies (frequency range: 50~540 kHz) were used for our PAT reconstruction.

For the experiments,^{12} pulsed light from a Nd: YAG laser (wavelength: 532nm, pulse duration: 3–6ns) were coupled into the phantom via an optical subsystem and generated acoustic signals. The transducer (1MHz central frequency) and phantom were immersed in a water tank. A rotary stage rotated the receiver relative to the center of the tank. The incident optical fluence was controlled below 10mJ/cm^{2} and the incident laser beam diameter was 5.0cm. For the first two experiments, we embedded two objects with a size ranging from 2.0–5.5 mm in diameter in a 50.8 or 40.0 mm-diameter solid cylindrical phantom. We then immersed the object-bearing solid phantom into a 110.6 mm-diameter water background. The phantom materials used consisted of Intralipid as scatterer and India ink as absorber with Agar powder (1–2%) for solidifying the Intralipid and India ink solution. The background phantom had *µ _{a}*=0.01 mm

^{-1}and µ

*′*=1.0mm

_{s}^{-1}while the two targets had

*µ*=0.03 mm

_{a}^{-1}and µ

*′*=2.0 mm

_{s}^{-1}for test 1, and

*µ*=0.07 mm

_{a}^{-1}and µ

*′*=3.0mm

_{s}^{-1}for test 2. In the next two experiments, we placed a single-target-containing phantom into the water, aiming to test the capability of resolving target having different optical contrasts relative to the background phantom. The target size was 1.0 and 2.0mm in diameter for tests 3 and 4, respectively. The target had

*µ*=0.03 mm

_{a}^{-1}and µ

*′*=2.0 mm

_{s}^{-1}for test 3, and

*µ*=0.015 mm

_{a}^{-1}and µ

*′*=2 mm

_{s}^{-1}for test 4. In the image reconstructions for the four experiments, we assumed scattering coefficient known as constant (1.0mm

^{-1}). The initial guesses of optical absorption coefficient for the target(s) and background medium were 0.02mm

^{-1}and 0.01mm

^{-1}, respectively. Although a single transducer is used, the transducer has a bandwidth that allows us to use multiple frequencies by simply Fourier transforming the detected time domain acoustic signals. In this work, 50 frequencies (frequency range: 50~540 kHz) were used for our PAT reconstruction. It required about 30 minutes to finish the two-step reconstruction computation.

## 3. Results and discussion

The results from simulated data are shown in Fig.1 where Fig. 1(a) provides the distribution of optical fluence, Fig. 1(b) presents the reconstructed absorbed energy density image using the existing PAT algorithm, and Fig. 1(c) displays the recovered absorption coefficient image with the regularized Newton method. We can see from Fig. 1(c) that absorption coefficient image can be recovered quantitatively. It is also observed from Figs. 1(a) and 1(b) that the influence of the inhomogeneous distribution of photon density on the PAT reconstruction is apparent. There is no linear relation existing between the absorbed energy density and optical absorption coefficient even if the incident distributed source is homogeneous, as demonstrated by Figs. 1(b) and 1(c).

The results from the first two sets of experiments are shown in Fig. 2 where Figs. 2(a) and 2(b) present the reconstructed absorption coefficient images of two objects having a size of 2.0 and 3.0mm (test 1), and 5.5mm (test 2) in diameter, respectively, while the recovered absorbed energy density maps for experiments 1 and 2 are also plotted in Figs. 2(c) and 2(d) for comparison. We see that the objects in each case are clearly detected. As shown in Table 1, the recovered absorption coefficients of the target and background are quantitative compared to the exact values for both experiments. By estimating the full width half maximum (FWHM) of the absorption coefficient profiles, the recovered object size was found to be 1.8, 2.7, and 5.0 mm, which is also in good agreement with the actual object size of 2.0, 3.0, and 5.5 mm for experiments 1 and 2. The reconstructed absorption coefficient images for experiments 3 and 4 are shown in Figs. 3(a) and 3(b). We immediately see that the different optical contrast levels of the objects relative to the background are quantitatively resolved.

It is important to note that our reconstruction method does not need any calibration procedure due to the use of relative incident laser source strength and normalized absorbed energy density distribution where the normalization is performed by simply dividing the absorbed energy density at each nodal location with the maximum absorbed energy density (i.e., Φ=Φr)(r)/_{Φ}). An optimization scheme was then applied to search for the boundary conditions coefficient, *α* and the relative source strength as described previously.^{15} As such, the reconstruction of optical properties with our algorithms does not depend on the absolute values of absorbed energy density and optical fluence as well as the boundary parameter. For example, even though the values/scales of the absorbed energy density for experiments 1 and 2 are very different as shown in Figs. 2(c) and 2(d), the algorithm is still able to recover the absorption coefficient images quantitatively in terms of the location, size, and absorption coefficient value of the objects. In addition, our method is able to resolve the issue of negative absorbed energy density values often seen in conventional PAT. For our previous methods^{12}, the negative values must be specified as zero, which may affect the quantitative accuracy of the recovered absorption coefficient images. In this study we demonstrate experimental evidence that it is possible to obtain absolute optical absorption coefficient image using photoacoustic tomography coupled with diffusion equation based regularized Newton method. The methods described are able to quantitatively reconstruct absorbing objects with different sizes and optical contrast levels.

## Acknowledgments

This research was supported in part by a grant from the National Institutes of Health (NIH) (R01 CA90533).

## References and links

**1. **R.G.M. Kollkman, W. Steenbergen, and T. G. van Leeuwen, “Reflection mode photoacoustic measurement of seed of sound,” Opt. Express **15**, 3291–3300 (2007). [CrossRef]

**2. **R. I. Siphanto, K. K. Thumma, R.G.M. Kolkman, T.G. van Leeuwen, F.F.M. de mul, J.W. van Neck, L.N.A. van Adrichem, and W. Steenbergen, “Serial noninvasive photoacoustic imaging of neovascularization in tumor angiogenesis,” Opt. Express **13**, 89–95 (2005). [CrossRef] [PubMed]

**3. **Z. Chen, Z. Tang, and W. Wan, “Photoacoustic tomography imaging based on a 4f acoustic lens imaging systems” Opt. Express **15**, 4966–4976 (2007). [CrossRef] [PubMed]

**4. **Z. Yuan, Q. Zhang, and H. Jiang, “Simultaneously reconstruction of acoustic and optical properties of heterogeneous medium by quantitative photoacoustic tomography,” Opt. Express **14**, 6749–6753 (2006). [CrossRef] [PubMed]

**5. **R.A. Kruger, D. Reinecke, and G. Kruger, “Thermoacoustic computed tomography-technical considerations”, Med. Phys. **26**, 1832–1837(1999). [CrossRef] [PubMed]

**6. **G. Paltauf, J. Viator, S. Prahl, and S. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging”, J. Acoust. Soc. Am. **112**, 1536–1544(2002). [CrossRef] [PubMed]

**7. **S.J. Norton and T. Vo-Dinh, “Optoacoustic diffraction tomography: analysis of algorihms”, J. Opt. Soc. Am. A **20**, 1859–1866(2003). [CrossRef]

**8. **A.A. Karabutov, E. Savateeva, and A. Oraevsky, “Imaging of layered structures in biological tissues with opto-acoustic front surface transducer”, Proc. SPIE **3601**, 284–295(1999). [CrossRef]

**9. **C.G.A. Hoelen, F.F. de Mul, R. Pongers, and A. Dekker, “Three-dimensional photoacoustic imaging of blood vessls in tissue”, Opt. Lett. **23**, 648–650 (1998). [CrossRef]

**10. **H. Jiang, Z. Yuan, and X. Gu, “Spatially varying optical and acoustic property reconstruction using finite element-based photoacoustic tomography,” J. Opt. Soc. Am. A **23**, 878–888 (2006). [CrossRef]

**11. **J. Ripoll and V. Ntziachristos, “Quantitative point source photoacoustic inversion formulas for scattering and absorbing medium,” Phys. Rev. E **71**, 031912 (2005). [CrossRef]

**12. **Z. Yuan and H. Jiang, “Quantitative photoacoustic tomography: recovery of optical absorption coefficient map of heterogeneous medium,” Appl. Phys. Lett. **88**, 231101 (2006). [CrossRef]

**13. **B. Cox, S. Arridge, K. Kostli, and P. Beard, “2D quantitative photoacoustic image reconstruction of absorption distributions in scattering medium using a simple iterative method,” Appl. Opt. **45**, 1866–1875 (2006). [CrossRef] [PubMed]

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

**15. **N. Iftimia and H. Jiang, “Quantitative optical image reconstruction of turbid media by use of direct-current measurements,” Appl. Opt. **39**, 5256–5261 (2000). [CrossRef]

**16. **P. Yalavarthy, H. Dehghani, B. Pogue, and K.D. Paulsen, “Wight-matrix structured regularization provide optimal generalized least-square in diffuse optical tomography,” Med. Phys. **34**, 2085–2098 (2007). [CrossRef] [PubMed]

**17. **Z. Yuan and H. Jiang, “Three-dimensional finite element-based photoacoustic tomography: Reconstruction algorithm and simulations,” Med. Phys. **34**, 538–546 (2007). [CrossRef] [PubMed]