Listening to light scattering in turbid media: quantitative optical scattering imaging using photoacoustic measurements with one-wavelength illumination

Conventional biomedical photoacoustic tomography (PAT) reconstruction methods are mostly qualitative based [1–3], which provide only the distribution of absorbed light energy density that is the product of the optical absorption coefficient and local optical fluence within the irradiated medium. Several methods recently reported suggest that it is possible to recover absorption coefficient maps from the PAT-derived absorbed energy density [4–12]. All these methods for quantitative PAT (qPAT) require a model-based reconstruction algorithm to firstly identify the map of absorbed energy density or initial pressure. Then under various assumptions, the absorption coefficient is deduced from the recovered absorbed energy density/initial pressure through photon diffusion equation or a measurement technique that can capture the distribution of exciting fluence [4–12]. However, the recovery of optical absorption coefficient based on the model-based method depends on the accuracy of the distribution of absorbed energy density and a parameter calibration procedure, whereas the experimental method uses invasive measurements of optical fluence, which is impractical for clinical applications. In particular, for most of the previous studies, the scattering coefficient was assumed a constant [13].

To achieve complete and true qPAT, the recovery of scattering coefficient must be addressed because distribution of scattering coefficient in biological tissue is typically heterogeneous and light scattering dominates over absorption by an order of magnitude or more [14]. Light scattering not only affects the accuracy of the recovered absorption coefficient if it is assumed a constant, but scattering coefficient itself is a very important intrinsic tissue parameter. The exploration of its recovery will open a new avenue to realize complete qPAT by providing spatially resolved cellular morphology such as cellular density and crowding for accurate diagnostic-decision making. While no successful experiments have been reported so far, previous attempts using simulated data showed that scattering coefficient could be recovered photoacoustically from multi-source illumination if the optical properties of background medium were assumed known [15]. The listening-to-light-scattering (LLS) method described here eliminates such assumptions and is able to provide quantitative recoveries of both scattering and absorption coefficients using photoacoustic measurements with one-wavelength illumination.

The developed LLS method is based on two steps. The first step is to directly extract both absorption coefficient and absorbed energy density distributions using PAT from the photoacoustic measurements along the boundary. Then the knowledge of the distributions of absorption coefficient and absorbed energy density is used to recover absorption and scattering coefficients simultaneously through the photon diffusion equation.

2.1. Methods

2.1.1. Reconstruction of absorbed energy density and absorption coefficient by PAT

As the first step, the absorbed energy density is recovered by a model-based PAT reconstruction algorithm. The core procedure of the algorithm can be described by the following equations responsible for the forward and inverse solutions, respectively [10, 16],

$$\begin{eqnarray}&&{{\nabla }^{2}}p(r,\omega )+k_{0}^{2}p(r,\omega )=i{{k}_{0}}\frac{{{c}_{0}}\beta E(r)}{{{C}_{p}}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&({{{\boldsymbol{J}} }^{T}}{\boldsymbol{J}} +\lambda {\boldsymbol{I}} )\Delta {\mathbf{\chi }}={{{\boldsymbol{J}} }^{T}}({{{\boldsymbol{p}} }^{o}}-{{{\boldsymbol{p}} }^{c}}),\end{eqnarray} \tag{ 2 }$$

in which p is the pressure wave; ${{k}_{0}}=\omega /{{c}_{0}}$ is the wave number described by the angular frequency, $\omega$ and the speed of acoustic wave in the medium, ${{c}_{0}}$; $\beta$ is the thermal expansion coefficient; C_p is the specific heat; E is absorbed energy density that is the product of absorption coefficient, ${{\mu }_{a}}$ and optical fluence, $\Phi$ (i.e., $E={{\mu }_{a}}\Phi$); ${{{\boldsymbol{p}} }^{o}}={{\left( p_{1}^{o},p_{2}^{o},...,p_{M}^{o} \right)}^{T}}$, ${{{\boldsymbol{p}} }^{c}}={{\left( p_{1}^{c},p_{2}^{c},...,p_{M}^{c} \right)}^{T}}$, where $p_{i}^{o}$ and $p_{i}^{c}$ are observed and computed complex acoustic field data for i = 1, 2,..., M boundary locations; $\Delta {\boldsymbol{ \chi }}$ is the update vector for E; J is the Jacobian matrix formed by ∂p/∂E at the boundary measurement sites; λ is a Levenberg–Marquardt regularization parameter and I is the identity matrix. The image formation task here is to update E via iterative solution of equations (1) and (2) so that a weighted sum of the least squared error between computed and measured acoustic data along the boundary domain can be minimized.

Similarly, ${{\mu }_{a}}$ can also be directly reconstructed through equation (1) based on the measured p along the boundary. In this case, the regularization based Newton method is used to update the initial guess of absorption coefficient ${{\mu }_{a,0}}$ based on the following equations (3) and (4)

$$\begin{eqnarray}&&({{\Im }^{T}}\Im +\lambda {\boldsymbol{I}} )\Delta {\boldsymbol{ \mu }} ={{\Im }^{T}}({{{\boldsymbol{p}} }^{o}}-{{{\boldsymbol{p}} }^{c}}),\quad {\boldsymbol{ \mu }} ={{{\boldsymbol{ \mu }} }_{a,0}}+\eta \Delta {\boldsymbol{ \mu ,}} \end{eqnarray} \tag{ 3 }$$

where $\Delta {\boldsymbol{ \mu }}$ is the update vector for ${{\mu }_{a}}$, $\Im$ is the Jacobian matrix formed by $\partial p/\partial {{\mu }_{a}}$ at the boundary measurement sites, and $\eta$ is calculated from the backtracking line search. In consideration of E = ${{\mu }_{a}}$ $\Phi$, the elements in Jacobian matrix $\Im$ are determined by,

$$\begin{eqnarray}\frac{\partial {{p}_{i}}}{\partial {{\mu }_{a,j}}}=\begin{array}{ccccccccccccccc} \underset{k=1}{\overset{N}{\mathop \sum }}\,\left( \frac{\partial {{p}_{i}}}{\partial {{E}_{k}}}\frac{\partial {{E}_{k}}}{\partial {{\mu }_{a,j}}} \right) & (i=1,2...M;j=1,2...N) \\ \end{array},\end{eqnarray} \tag{ 4 }$$

in which N is the total number of unknown parameters. The derivative $\partial E/\partial {{\mu }_{a}}$ in equation (4) is further written,

$$\begin{eqnarray}\frac{\partial {{E}_{k}}}{\partial {{\mu }_{a,j}}}=\left( \begin{array}{ccccccccccccccc} {{\Phi }_{j}}+{{\mu }_{a,j}}(\partial {{\Phi }_{j}}/\partial {{\mu }_{a,j}}) & if\begin{array}{ccccccccccccccc} \begin{array}{ccccccccccccccc} {} & (k=j) \\ \end{array} \\ \end{array} \\ {{\mu }_{a,k}}(\partial {{\Phi }_{k}}/\partial {{\mu }_{a,j}}) & if\begin{array}{ccccccccccccccc} {} & (k\ne j) \\ \end{array} \\ \end{array} \right.,\end{eqnarray} \tag{ 5 }$$

The sensitivities of $\partial p/\partial E$ and $\partial \Phi /\partial {{\mu }_{a}}$ in each iteration can be calculated using adjoint method from both equation (1) and the photon diffusion equation (see equation (6) below), respectively.

2.1.2. Reconstruction of absorption and scattering coefficients by LLS

In the second step, we use E recovered by the first step as the measured field data, and implement an iterative procedure for reconstructing the heterogeneous distributions of both the absorption and scattering coefficients starting from a spatially-varying initial guess of absorption coefficient (obtained by the first step) and a homogeneous initial guess of scattering coefficient, based on the following photon diffusion equation associated with the Robin boundary conditions [14, 17–19],

$$\begin{eqnarray}&&\nabla \cdot D(r)\nabla \Phi (r)-{{\mu }_{a}}(r)\Phi (r)=-S(r),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&-D\nabla \Phi \cdot n=\alpha \Phi ,\end{eqnarray} \tag{ 7 }$$

in which D(r) is the diffusion coefficient, $\alpha$ is a coefficient related to the internal reflection at the boundary, and S(r) is the normalized source term. The diffusion coefficient can be written as $D=1/(3({{\mu }_{a}}+{{\mu }_{s}}\prime )),$ where ${{\mu }_{s}}\prime$ is the reduced scattering coefficient. Here the goal is to recover the optical properties using the measured E inside the whole problem domain. In this case, the objective function for the minimization statement is given,

$$\begin{eqnarray}&&F=\underset{i=1}{\overset{N}{\mathop \sum }}\,{{\left( E_{i}^{m}-E_{i}^{c} \right)}^{2}}.\end{eqnarray} \tag{ 8 }$$

The objective function can be minimized by specifying Min: F = 0. This is a typical optimization problem, in which we are particularly interested in finding the optical property $\zeta$ that makes Min: F close to zero. Following the Taylor series expansion method, we obtain the approximated optical property $\zeta$ (absorption and scattering coefficients) from nearby point ${{\zeta }_{0}}$ ($\Delta \zeta =\zeta -{{\zeta }_{0}}$),

$$\begin{eqnarray}&&\frac{\partial F}{\partial \zeta }=\frac{\partial F}{\partial \zeta }({{\zeta }_{0}})+\Delta \zeta \frac{{\rm{d}}}{{\rm{d}}\zeta }\left[ \frac{\partial F}{\partial \zeta }({{\zeta }_{0}}) \right]+...\end{eqnarray} \tag{ 9 }$$

If the higher order terms are ignored and $\frac{\partial F}{\partial \zeta }=0$ assumed, equation (9) is rewritten,

$$\begin{eqnarray}&&\zeta ={{\zeta }_{0}}-{{\left\{ \frac{{\rm{d}}}{{\rm{d}}\zeta }\left[ \frac{\partial F}{\partial \zeta }({{\zeta }_{0}}) \right] \right\}}^{-1}}\frac{\partial F}{\partial \zeta }({{\zeta }_{0}}).\end{eqnarray} \tag{ 10 }$$

The iterative format for equation (10) can be specified as

$$\begin{eqnarray}&&{{\zeta }_{i}}={{\zeta }_{i-1}}-{{\left\{ \frac{{\rm{d}}}{{\rm{d}}\zeta }\left[ \frac{\partial F}{\partial \zeta }({{\zeta }_{i-1}}) \right] \right\}}^{-1}}\frac{\partial F}{\partial \zeta }({{\zeta }_{i-1}}).\end{eqnarray} \tag{ 11 }$$

Based on equation (8), we can solve for the first- and second-order derivatives of F as,

$$\begin{eqnarray}&&\frac{\partial F}{\partial \zeta }=2{{\left( \frac{\partial {{E}^{c}}}{\partial \zeta } \right)}^{T}}({{E}^{c}}-{{E}^{m}}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\frac{{{\partial }^{2}}F}{\partial {{\zeta }^{2}}}=2{{\left( \frac{\partial {{E}^{c}}}{\partial \zeta } \right)}^{T}}\frac{\partial {{E}^{c}}}{\partial \zeta }+2{{\left( \frac{{{\partial }^{2}}{{E}^{c}}}{\partial {{\zeta }^{2}}} \right)}^{T}}({{E}^{c}}-{{E}^{m}}).\end{eqnarray} \tag{ 13 }$$

The contribution from the higher order derivative term in equation (13) is small and often discarded. Then inserting equation (13) into equation (11), we get

$$\begin{eqnarray}&&{{\zeta }_{i}}={{\zeta }_{i-1}}+{{\left\{ 2{{\left( \frac{\partial {{E}^{c}}}{\partial \zeta } \right)}^{T}}\frac{\partial {{E}^{c}}}{\partial \zeta } \right\}}^{-1}}2{{\left( \frac{\partial {{E}^{c}}}{\partial \zeta } \right)}^{T}}({{E}^{m}}-{{E}^{c}}),\end{eqnarray} \tag{ 14 }$$

in which ∂E/∂ζ is the Jacobian matrix G generated inside the whole problem domain. It should be noted that the Hessian matrix ${{{\mathbf{G}}}^{T}}{\mathbf{G}}$ in equation (14) is always ill-conditioned, which makes the iteration process unstable. A typical way to stabilize the inversion problem is through regularization to make ${{{\mathbf{G}}}^{T}}{\mathbf{G}}$ more diagonally dominant. So the ultimate iterative equation for updating the optical properties becomes

$$\begin{eqnarray}&&({{{\boldsymbol{G}} }^{T}}{\boldsymbol{G}} +\lambda \prime {\boldsymbol{I}} )\Delta {\boldsymbol{ \zeta }} ={{{\boldsymbol{G}} }^{T}}({{{\boldsymbol{E}} }^{m}}-{{{\boldsymbol{E}} }^{c}}),\end{eqnarray} \tag{ 15 }$$

in which $\lambda \prime$ is a scalar and $\Delta \zeta ={{(\Delta {{D}_{1}},\Delta {{D}_{2}},...,\Delta {{D}_{N}},\Delta {{\mu }_{a,1}},\Delta {{\mu }_{a,2}},...,\Delta {{\mu }_{a,N}})}^{T}}$. The sub-matrix in G can be calculated from equation (5) for ${{\mu }_{a}}$ and equation (16) for D, respectively,

$$\begin{eqnarray}&&\frac{\partial {{E}_{k}}}{\partial {{D}_{j}}}={{\mu }_{a,k}}(\partial {{\Phi }_{k}}/\partial {{D}_{j}}).\end{eqnarray} \tag{ 16 }$$

Thus the absorption and scattering coefficients are reconstructed through the iterative procedures described by equations (6) and (15).

2.2. Simulation and experimental tests

The image formation process described above is first tested using simulated data. The test geometry is shown in figure 1(a) where a two-dimensional circular background region (50.0 mm in diameter) contained two circular targets (5.0 mm in diameter each). The optical properties for the bottom left and top right targets were ${{\mu }_{a}}$ = 0.1 mm⁻¹/${{\mu }_{s}}\prime$ = 1.0 mm⁻¹ (scattering coefficient) and ${{\mu }_{a}}$ = 0.01 mm⁻¹/${{\mu }_{s}}\prime$ = 2.0 mm⁻¹, respectively, while the optical properties for the background medium were ${{\mu }_{a}}$ = 0.01 mm⁻¹/${{\mu }_{s}}\prime$ = 1.0 mm⁻¹. For this simulation test, a uniform distributed area source was utilized to illuminate the whole imaging volume from the top surface, as plotted in figure 1. The simulation results are shown in figure 1 where figures 1(b)–(d) present the exact distribution of energy density E (b), and the ${{\mu }_{a}}$ (c) and ${{\mu }_{s}}\prime$ (d) images reconstructed by the LLS method. The images shown in figures 1(c) and (d) and the property profiles presented in figures 1(e) and (f) indicate that both the absorption and scattering images can be simultaneously recovered and quantitatively separated with zero cross-talk between the two parameters. It is also observed from figure 1(b) that the influence of the inhomogeneous distribution of optical fluence on the absorbed energy density is apparent, where low values of energy density appear in the center area of the target. However, figures 1(c) and (d) indicate that the LLS method has the capability to resolve the ${{\mu }_{a}}$ and ${{\mu }_{s}}\prime$ maps from the heterogeneous distribution of local E completely.

Zoom In Zoom Out Reset image size

We then present a series of carefully designed phantom experiments to provide solid evidence that both the absorption and scattering images of heterogeneous scattering media can be concurrently reconstructed by the LLS method when the absorption coefficient distribution obtained by conventional PAT is incorporated into the LLS reconstruction procedure.

For the experimental setup shown in figure 2, a pulsed light beam from a Nd:YAG laser (wavelength: 532 nm, pulse duration: 3–6 ns) was sent to the phantom via an optical subsystem and generated acoustic signals. A transducer (1 MHz central frequency) and the 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 10 mJ cm⁻² and the incident laser beam diameter was 5.0 cm at the phantom surface. For the first two experiments, we embedded three objects with 4.0 mm in diameter in a solid cylindrical phantom, in which each object had different absorption and scattering contrasts, to validate if both the parameters can be recovered simultaneously and the cross-talk between the two parameters can be resolved. We then immersed the object-bearing solid phantom into the 80.0 mm-diameter water tank. 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 ${{\mu }_{a}}$ = 0.01 mm⁻¹ and ${{\mu }_{s}}\prime$ = 1.0 mm⁻¹ while the top left target had ${{\mu }_{a}}$ = 0.02 mm⁻¹ and ${{\mu }_{s}}\prime$ = 1.0 mm⁻¹, the top right target had ${{\mu }_{a}}$ = 0.02 mm⁻¹ and ${{\mu }_{s}}\prime$ = 2.0 mm⁻¹, and the bottom right target had ${{\mu }_{a}}$ = 0.01 mm⁻¹ and ${{\mu }_{s}}\prime$ = 2.0 mm⁻¹ for test 1. To verify the capability of resolving targets having different optical contrasts relative to the background phantom, the following optical properties were used for test 2: ${{\mu }_{a}}$ = 0.04 mm⁻¹ and ${{\mu }_{s}}\prime$ = 4.0 mm⁻¹ for the top left target, ${{\mu }_{a}}$ = 0.01 mm⁻¹ and ${{\mu }_{s}}\prime$ = 4.0 mm⁻¹ for the top right target, and ${{\mu }_{a}}$ = 0.04 mm⁻¹ and ${{\mu }_{s}}\prime$ = 1.0 mm⁻¹ for the bottom middle target.

Zoom In Zoom Out Reset image size

The reconstructed results from the first two sets of experiments are shown in figures 3 and 4, where figures 3(b) and 4(b) present the reconstructed E images for tests 1 and 2, respectively, while the recovered ${{\mu }_{a}}$ maps using PAT for experiments 1 and 2 are plotted in figures 3(c) and 4(c). The recovered ${{\mu }_{a}}$ and ${{\mu }_{s}}\prime$ images for the two experiments by LLS method are presented in figures 3(d) and 4(d) and figures 3(e) and 4(e), respectively. The ${{\mu }_{a}}$ and ${{\mu }_{s}}\prime$ images reconstructed by LLS for the hair-containing phantom are plotted in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results presented in figures 3 and 4 clearly suggest that simultaneous reconstruction of both absorption and scattering images of heterogeneous scattering media can be obtained using LLS. In particular, we see that the target with absorption contrast appears only in the absorption images and the target with scattering contrast shows only in the scattering maps, giving zero cross-talk between the two parameters. We also found that when the targets had both absorption and scattering contrasts, the images show better quality in terms of the recovered size and quantitative values of optical properties of the targets. However, when the target was a pure absorber or a pure scatter, the recovered images show underestimated optical properties, which is clearly illustrated by the quantitative property profiles shown in figures 3(f)–(i). Additionally, the scattering coefficient profile given in figure 3(i) indicates that the size of the targets with pure scattering contrast only was significantly underestimated. Taking a close examination of the absorption images recovered by LLS and PAT (figures 3 and 4), we observe that the absorption images recovered by the two methods give similar quality. This suggests that we could use the absorption images recovered by PAT or LLS in practice, while the scattering images could only be recovered by the LLS method. Figure 5 shows that the ${{\mu }_{a}}$ and ${{\mu }_{s}}\prime$ images of human hair phantom are clearly recovered by LLS with ultrasound resolution, which is impossible for pure optical imaging methods such as diffuse optical tomography (DOT).

In addition, it is noted from the absorption image in figure 5 that structural broken phenomenon shows up for this human hair phantom test. We think this is due to the inhomogeneous light illumination during the experiment though we try to generate a uniform distribution of laser sources, which will definitely affect the sensitivity (SNR) of the photoacoustic imaging systems. As a result, the imaging resolution for optical absorption coefficient can be degraded for part of the human hair phantom. For example, both the recovered size and the absorption coefficient for part of the hair can be changed.

Further it should be pointed out that the existence of a solution to an inversion problem depends on not only just the uniqueness, but also other issues including the quality of the data (SNR) and the statistical behavior of the problems (complex in vivo data; having only absorption/scattering contrast or both) [20]. For example, it is observed from figures 3 and 4 that when the heterogeneities/objects have contrast in both the absorption and scattering properties, the optical properties can be recovered with highly quantitative accuracy, which will not be strongly affected by the noise. However, when targets only have absorption or scattering contrast, the image quality will be degraded. This is largely due to the fact that unlike the absorption coefficient which is strongly related to the absorbed energy density distribution, the scattering coefficient affects the absorbed energy density indirectly through optical fluence. In particular, for pure absorption cases, the distribution of the optical fluence will also have some effect on the reconstruction accuracy of optical absorption coefficient.

Based on the developed reconstruction method, we have repeatedly shown using controlled phantom experiments that absorption and scattering coefficient images can be simultaneously recovered when conventional Tikhonov regularization is combined with a normalization scheme. However, for in vivo tests using measurements with strong noise, the scheme may fail to recover the absorption and scattering coefficients for some extreme or complex cases.

In future, for very noisy in vivo test, we plan to use multispectral PAT reconstruction methods, which are able to allow us to recover both scattering and absorption-related parameters with enhanced imaging quality [21]. General multispectral PAT methods are needed in order to minimize the reconstruction errors with increased photoacoustic measurements.

In summary, we have presented solid experimental evidence that quantitative absorption and scattering images of heterogeneous scattering media can be recovered simultaneously with ultrasound resolution by use of photoacoustic data. In particular, the results from controlled phantom experiments have shown that the interparameter cross-talk often seen in DOT is eliminated using the LLS scheme. Further evaluation of LLS using in vivo clinical data is underway in our laboratory and the results will be reported elsewhere.

This research was supported by SRG2013-00035-FHS Grant and MYRG2014-00093-FHS Grant from University of Macau in Macau.