1.

Introduction

Photoacoustic (PA) imaging has been intensively studied recently because of its promise for high resolution and intrinsic optical contrast.¹ Efforts have been directed to quantitative estimation of optical properties (absorption coefficient ${\mu }_a$ and scattering coefficient ${\mu }_s$), and the Grüneisen parameter distributions using photoacoustic tomography (PAT), which may significantly enhance clinical and biomedical applications of PA imaging. This is usually termed quantitative photoacoustic tomography (qPAT), which involves two steps. The first step is a well-studied acoustic inverse problem, which aims to reconstruct PA initial pressure distributions using recorded acoustic data generated by a wide field-pulsed laser.^2–4 The second step is to estimate optical properties ${\mu }_a$ and ${\mu }_s$. This is essentially a rather challenging optical inverse problem. Quite a number of approaches have been proposed for the second step of qPAT. Studies were first focused on ${\mu }_a$ estimation,^5–7 and later extended to include scattering.^8–10 Multiple-wavelength^11,12 or optical sources^13–16 were also proposed to mitigate the absorption-scattering nonuniqueness problem. Diffuse optical measurements were also used to estimate ${\mu }_a$ in PAT.^17,18 A thorough review on qPAT is provided by Cox et al.¹⁹

The fixed-point iterative scheme for estimating optical absorption distribution was first proposed by Cox et al.⁶ With this method, absorption maps are estimated in each iteration and then used to update fluence distributions for the next iteration. With simulations based on a diffusion-based finite-element model, the authors demonstrated that the algorithm converges fast and the absorption distribution could be reconstructed accurately. Jetzfellner et al.²⁰ examined the performance of this scheme with experimental data. With a tissue-mimicking phantom, PA imaging was conducted using circular illumination with a high intensity near-infrared-pulsed laser. The authors argued the iterative method is sensitive to background optical properties and diverges over iterations.

Recently, Harrison et al.²¹ extended the aforementioned iterative method⁶ to a least-squares fixed-point iterative method for reconstructing absorption maps in multiple-illumination photoacoustic tomography (MIPAT). Instead of using single-optical source, a number of optical illuminations evenly distributed around the subject in the imaging plane are employed in MI-PAT (as is shown in Fig. 1). With the proposed method, convergence of the inversion is significantly improved when multiple illuminations are utilized. This approach does not require inversion of large Jacobian matrices, and it is efficient and stable. However, the aforementioned methods are restricted to the reconstruction of only absorption distributions.

Fig. 1

Schematic of the multiple-illumination photoacoustic tomography (MI-PAT).

In the present paper, we intend to extend our previous work to recover both absorption and scattering perturbation distributions in a known turbid media, which holds a typical value of background scattering in soft tissue. We accomplish this task in an iterative manner. For each iteration, absorption maps are first estimated with the least-squares fixed-point iterative method by Harrison et al.²¹ The optical fluence distribution is then updated, whereupon the scattering coefficient perturbation distribution is estimated. This procedure is then repeated multiple times until an acceptable error is reached.

2.

Method

2.3.

Reconstruction of the Diffusion Coefficient Perturbation Distribution

The relation between the measured data $p_i(\mathbf{r})$ and $[{\mu }_a(\mathbf{r}),D(\mathbf{r})]$ is nonlinear, because the optical fluence at location $\mathbf{r}$ is also a function of local optical properties. However, if we have an estimation of the absorption distribution $\widehat{{\mu }_a}(\mathbf{r})$, the problem is simplified as the contribution of absorption to $\widehat{{\mathrm{\Phi }}_i}(\mathbf{r})$ is known. Equation (2) is rewritten as

Eq. (4)

$$p_i(\mathbf{r})=\mathrm{\Gamma }(\mathbf{r})\widehat{{\mu }_a}(\mathbf{r})\widehat{{\mathrm{\Phi }}_i}[\widehat{{\mu }_a}(\mathbf{r}),D(\mathbf{r})].$$

To linearize the problem, we further decompose the local optical fluence at position $\mathbf{r}$ as $\widehat{{\mathrm{\Phi }}_i}(\mathbf{r})=\widehat{{\mathrm{\Phi }}_{i,o}}(\mathbf{r})+\delta \widehat{{\mathrm{\Phi }}_{i,a}}(\mathbf{r})+\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})$, which is a sum of known homogeneous background fluence ($\widehat{{\mathrm{\Phi }}_{i,o}}(\mathbf{r})$) due to known background absorption and scattering properties, and fluence perturbations due to absorption and scattering: $\delta \widehat{{\mathrm{\Phi }}_{i,a}}(\mathbf{r})$ and $\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})$.

Equation (4) now becomes

Eq. (5)

$$p_i(\mathbf{r})=\mathrm{\Gamma }(\mathbf{r})\widehat{{\mu }_a}(\mathbf{r})[\widehat{{\mathrm{\Phi }}_{i,o}}(\mathbf{r})+\delta \widehat{{\mathrm{\Phi }}_{i,a}}(\mathbf{r})+\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})].$$

In the above equation, the only unknown term is $\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})$, the fluence perturbation due to scattering diffusion coefficient abnormalities. In this study, an extrapolated boundary condition is used (with zero inward flux).²² We assume that $\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})$ is linearly related to $\delta D(\mathbf{r})$, which limits our study to the linear case such that $\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})\ll {\mathrm{\Phi }}_{i,o}(\mathbf{r})$, and ${\mathrm{\Phi }}_i\cong {\mathrm{\Phi }}_{i,o}$, then we have the following relation based on the first-order Born approximation:¹⁴

Eq. (6)

$$\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})=\int \frac{\delta D({\mathbf{r}}^{\prime })}{D_o}\nabla G_0(\mathbf{r},{\mathbf{r}}^{\prime })·\nabla {\mathrm{\Phi }}_{i,o}({\mathbf{r}}^{\prime })\mathrm{d}{\mathbf{r}}^{\prime }.$$

Here, $G_0(\mathbf{r},{\mathbf{r}}^{\prime })$ is Green’s function representing propagation from $\mathbf{r}$ to ${\mathbf{r}}^{\prime }$. $\delta {\mathrm{\Phi }}_{i,s}(r)$ can be estimated from Eq. (5):

Eq. (7)

$$\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})=\frac{p_k(\mathbf{r})}{\mathrm{\Gamma }\widehat{{\mu }_a}(\mathbf{r})}-\widehat{{\mathrm{\Phi }}_{i,o}}(\mathbf{r})-\delta \widehat{{\mathrm{\Phi }}_{i,a}}(\mathbf{r}).$$

The integral in Eq. (6) can be discretized as

Eq. (8)

$$\delta {\mathrm{\Phi }}_{i,s}(\mathbf{r})={\mathrm{\Sigma }}_n{W}_{\{ij\}n}^s\delta D({\mathbf{r}}_n),$$

where

Eq. (9)

$$W_{\{ij\}n}^s=\nabla G_0(\mathbf{r},{\mathbf{r}}_n^{\prime })·\nabla \mathrm{\Phi }({\mathbf{r}}_n^{\prime })·V/D_0.$$

$\nabla V$ is a volume element. $i$ ($i=\mathrm{1,2},\dots S$) is the index of optical sources, whereas $j$ is the index of optical fluence locations, with $j=\mathrm{1,2},\dots J$. The $\{ij\}$ elements are grouped together in this way because they could be rasterized into a single index $k_{ij}=i+jS$.

We write Eq. (7) in a matrix form as

Eq. (10)

$${\mathrm{\Phi }}_{\mathbf{sc}}=\mathbf{W}\mathbf{u},$$

where $\mathbf{W}$ is the Jacobian matrix ($SJ\times N$), or sensitivity matrix, and $\mathbf{u}$ is the $N\times 1$ vector which contains the unknown diffusion coefficient perturbation map $\delta D(\mathbf{r})$ to estimate. This is a linearized inverse problem and can be solved with a number of techniques. We use the least-squares minimization with total variation regularization,

Eq. (11)

$$\mathbf{u}={\arg \min }_{\mathbf{u}}{\Vert {\mathrm{\Phi }}_{\mathbf{sc}}-\mathbf{W}\hat{\mathbf{u}}\Vert }_2^2+\lambda {|\hat{\mathbf{u}}|}_{\mathrm{TV}}.$$

$\lambda$ in the above equation serves as the regularization parameter for inversion. To determine the value of $\lambda$, we use the classic $L$-curve technique to balance the resulting errors in the first and second term.²³

We propose to repeat this procedure with multiple iterations to accomplish the imaging task until a preset tolerance is reached following the flowchart in Fig. 2. Before entering the iteration loop, we set the initial guess of the absorption as zero. In each iteration, we first update the fluence distribution with the estimated $[{\mu }_a,D]$ from the last iteration. Then, we reconstruct $D$ with the updated flunce distribution, whereupon a new version of optical fluence is generated, based on which absorption is again modified if a tolerance is not yet reached. The iterative reconstruction is repeated until a satisfactory result is reached. There are two types of iterations involved in this algorithm. To distinguish them, we call the repetition of consecutively reconstructing $[{\mu }_a^{(n)},D^{(n)}]$ an iteration, and the repetitive steps for estimating ${\mu }_a^{(i)}$ the ${\mu }_a$-iteration.

Fig. 2

Flow chart of the iterative method.

3.

Simulation

Simulation studies were conducted to examine the performance of the proposed scheme. All simulations were conducted on a region of interest composed of a two-dimensional (2-D) $40\times 40$ mesh grid system, whose dimension is set to be $20\mathrm{mm}\times 20\mathrm{mm}$ unless otherwise stated. To mitigate inversion crime, noise was added and the synthetic data were generated on a larger grid system ($43\times 43$) and down-sampled to $40\times 40$ with the cubic interpolation for inversion. 16 optical sources were positioned 3-mm backward from the object to validate the diffusion light propagation model. The background optical properties of the object were taken as the physiologically realistic values in biological soft tissue. Absorption and the reduced scattering coefficient are 0.1 and $100{\mathrm{cm}}^{-1}$, respectively.

We first tested our algorithm with regular overlapping features, as is shown in Fig. 3. Two rectangular features with absorption perturbations (10% on the left and 5% on the right) and diffusion coefficient anomalies (5% on the left and 10% on the right) were positioned in the center of the field-of-view. White noise with a normal distribution corresponding to an average signal-to-noise (SNR) of $\sim 30\mathrm{dB}$ was added to “reconstructed” PAT images. Figure 3(a) shows the PAT images due to two optical sources located on the top and bottom of the phantom, respectively. For inversion, we use zero as an initial guess for both absorption map and diffusion coefficient distribution. In each iteration, the ${\mu }_a$ iteration was conducted for 20 times. Reconstructed results with iteration 1, 2, 10, and 50 are employed to show improvement of estimation as more iterations are involved, as is shown in Fig. 3(b).

Fig. 3

Reconstruction of optical property distributions with the proposed method using 16 optical sources. (a) Synthetic photoacoustic data (initial pressure distribution) due to two optical sources located on top and bottom of the phantom, respectively. (b) Simulation results with iteration #0 (initial guess), #1, #2, #10, and #50 are shown to demonstrate that reconstruction results are improved with more iterations. Initial values were chosen to be zero for both absorption and diffusion coefficient distributions.

To examine the capability of the method to recover arbitrary features rather than regular shapes, features with both discrete strips with sharp edges and smoothed features were tested, as is shown in Fig. 4. In Fig. 4(a), sharp strip features with 10% absorption and scattering perturbations were reconstructed. We also tested four smoothed absorption perturbations (0.1, 0.09, 0.07, and $0.05{\mathrm{cm}}^{-1}$, respectively) and two diffusion coefficient anomalies (0.0083 and $0.0042{\mathrm{cm}}^{-1}$) with different spatial locations in the background. Synthetic data were generated with an SNR of 30 dB and 20 iterations were used for inversion. Values of the regularization parameter were different for the example to honor the feature characteristics.

Fig. 4

Reconstruction of sharp (a) and smoothed features (b) with the proposed method after 20 iterations. In (a), anomalies of both absorption and diffusion coefficient was set at 10% of the background value. Perturbations of 10%, 9%, 7%, and 5% of the background absorption, 20% and 10% of the background diffusion coefficient at various locations were used in (b). SNR of synthetic data was 30 dB.

Cross-talk between reconstructed maps of the two parameters is a concern. Artifacts due to diffusion coefficient features may appear in the reconstructed absorption map. Despite the fact that the diffusion coefficient depends on both absorption and scattering: $D(\mathbf{r})=1/3[{\mu }_a(\mathbf{r})+{\mu }_s^{\prime }(\mathbf{r})]$ locally, the scattering coefficient dominates the apparent $D$ values. If $D$ is not correctly reconstructed, the resulting fluence estimation may “contaminate” the reconstructed absorption map. We tried estimating two nonoverlapping features to examine cross-talk between the two properties in the recovered results (Fig. 5) in cases of single and multipleilluminations. Interestingly, whereas no artifacts due to absorption appeared in the reconstructed $D$, we did observe artifacts in the reconstructed absorption map where diffusion perturbations lie in both cases. However, the artifacts are alleviated within several iterations and were finally removed from the results when 16 optical sources were used. This might be because in each iteration, we update the absorption distribution, thus the dependence of $D$ on ${\mu }_a$ is significantly mitigated and is already negligible before reconstructing the scattering feature. On the other hand, while we estimate absorption maps, the contribution of scattering features brings about artifacts in the recovered ${\mu }_a$. But with more iterations, the optical fluence is updated repeatedly and approaches its true value, therefore, artifacts are removed. It is not surprising that with a single source, a unique solution for $[{\mu }_a,D]$ estimation is not achieved. This is due to the nonuniqueness problem.

Fig. 5

Reconstruction with nonoverlapping features to demonstrate capabilities of the algorithm to alleviate cross-talk between absorption and diffusion coefficient. Reconstruction results when using 1 and 16 illuminations are shown. Artifacts appeared in the reconstructed absorption map where diffusion coefficient anomalies locate. But with more iterations involved, the artifacts were removed when multiple-sources were used.

Figure 6 shows the reconstruction of a synthetic phantom with features that shares similarities with blood vessel vasculature using 16 sources. However, we would like to note that the primary purpose of this example is to challenge the inverse problem. The software phantom we used here may not be truly physiological, as certain wavelength bands optical absorption of blood may be too high to use a diffusion approximation. From left to right are the true model, reconstruction results with SNR of infinity (0 noise), $\sim 30$ and $\sim 20\mathrm{dB}$, respectively. We also show the cross-sectional profiles of the corresponding results in the center of the distributions in Fig. 7. With data that has an SNR of $\sim 30\mathrm{dB}$, absorption features are faithfully estimated. More artifacts appeared in the results computed with 20-dB data.

Fig. 6

Reconstruction of a synthetic phantom with features that shares similarities with blood vessel vasculature using 16 sources. The left column is the true $[{\mu }_a,D]$ model. Second to the last column are reconstruction results with average SNR of infinity (no noise), $\sim 30$ and $\sim 20\mathrm{dB}$.

Fig. 7

Cross-sectional profiles of the reconstruction results along the horizontal direction in the center of the true model in Fig. 6. All features of ${\mu }_a$ were faithfully reconstructed with the $\sim 30\text{-}\mathrm{dB}$ data. Artifacts appeared in the reconstruction with data that has an SNR of $\sim 20\mathrm{dB}$. Despite some artifacts, $D$ features were recovered with a $\sim 30\text{-}\mathrm{dB}$ SNR. Amplitude errors were unacceptable for data with $\sim 20\text{-}\mathrm{dB}$ SNR.

Figure 8 depicts the relation between reconstruction errors from the last example and the number of optical sources involved. Errors were calculated with the following equation: ²⁰

Fig. 8

Reconstruction errors with different numbers of optical sources.

Fast convergence was seen as a remarkable advantage of the least-squares fixed-point MIPAT approach.²¹ We plot the relative errors of reconstruction in Fig. 9 as an illustration of the overall convergence ability of the proposed scheme. Reconstruction errors decrease quickly and stabilize within the first 10 iterations.

Fig. 9

Normalized reconstruction errors with 50 iterations. Both ${\mu }_a$ and $D$ converge in a few iterations and stay stable after a large number of iterations.

4.

Conclusion and Discussion

We have proposed a consecutive reconstruction scheme to estimate optical absorption and diffusion coefficient perturbation distributions in a known turbid medium for MIPAT. With this approach, the absorption map is first estimated, whereupon the optical fluence is updated and then the diffusion coefficient distribution is reconstructed. The absorption distribution is recovered with a least-squares fixed-point iterative method. Numerical simulations demonstrated that both optical properties can be faithfully reconstructed with the presence of noise at an SNR level of $\sim 30\mathrm{dB}$. This method converges within only a number of iterations and is robust to cross-talk of overlapping and nonoverlapping features.

Conventional iterative methods for qPAT suffer from computational complexity, instability, and poor convergence performance. The least-squares fixed-point iterative methods²¹ to recover absorption in our previous report does not require inversion of large scale Hessian or Jacobian matrices, therefore, it is computationally efficient and stable. The reconstruction can be done on a resolution scale equivalent to the fine resolution obtainable by backprojection or other model-based inversion approaches. Convergence is substantially improved. These advantages enable the faithfully estimated absorption to serve as a compulsory guide for each iteration in the current method. Updated fluence estimate based on the reconstructed absorption map assisted accurate estimation of the scattering features. We did not study convergence of the algorithm with a much greater scattering background as described in Ref. 20 as we restrict our study to the scope of realistic biological soft tissue. Although various approaches have been proposed for quantitative reconstructions of optical distributions, most previous approaches were not sufficiently tolerant to realistic noise levels. Prior to experimental work, we argue it is essential to develop imaging strategies and algorithms for sufficient noise-robustness. In this manuscript, we demonstrate the ability to reconstruct both absorption and scattering distributions with realistic noise levels.

Our study has limitations. The result in this paper is based on the approximate 2-D light propagation model in the diffuse regime. Although many studies used the diffuse approximation for qPAT, the radiative transfer equation ²⁴ provides a better estimation. It is worth trying to extend our current study to more accurate theoretical models. Similar to most of the literature on this topic, we also assume ideal reconstruction of the PA images from measured acoustic data. Although one option to avoid this is to use the ultrasonic channel data as proposed by our previous report,¹⁶ one can also follow the framework by Saratoon et al.²⁵ to include acoustic reconstruction to study its influence on the final reconstruction results. Reconstruction of the scattering distributions still requires matrix inversion, as do all other diffuse optical tomographic (DOT) methods thus far. The matrix condition number is comparable with previous DOT/PAT approaches and the reconstruction can be done on a courser scale to accommodate regularization (data not shown). Resolution of scattering distributions may need to be courser than for estimated absorption distributions, which can be reconstructed with a resolution as fine as the transducers and algorithms will permit. We utilized this method as a preliminary example. We partially committed the “inverse crime” in this paper, which refers to the act to generate as well as to invert synthetic data with the same theoretical model or discretization. Since we use the same theoretical model for both the forward and inverse problems, we (1) generated synthetic data on a larger grid system (with an odd-number of grid) and then down-sampled the dataset with a smaller system (even number of grid points) for inversion; (2) added white noise with a normal distribution to all experimental data to mitigate this issue. Extension of our simulation work to experiments is planned as future work. As a preliminary work, recently we have reported a reflection-mode PA technique that can be used for optical sensing.²⁶ Our method requires inversion of the Jacobian matrix for estimation of the diffuse coefficient. Demand for computation resources increase significantly with Jacobian-matrix-based methods, including both memory and operation count,²⁵ therefore computation expense is a concern for the present framework. For example, despite the quick convergence of our method within only a few iterations, the computation time taken for each iteration with an $80\times 80$ grid system is $\sim 40\min$ with MATLAB 2010b (Mathworks Inc.) on a 64-bit 2.5-GHz Intel i5 PC system.

The main contribution of this paper is the demonstration of the potential to recover not only absorption but also scattering distributions with MIPAT, and to provide an algorithm with promising stability in moderately challenging noise. The study is meant to motivate experimental realization, which, however, is yet to be determined and optimized. Future work should aim to develop experimental embodiments of MIPAT and may require calibrating detected signals to known absorbed energy distributions, accounting for three-dimensional rather than 2-D light and sound propagations, acoustic attenuation, and transducer spatio-temporal impulse responses, among other factors.

If future experimental embodiments are successful, our approach could provide reconstruction of both absorption and scattering distributions in living subjects. Although DOT can provide such information, it does so at course resolution scales with poor depth-to-resolution ratios. MIPAT has the potential to quantitatively reconstruct absorption distributions with ultrasonic spatial resolution. The recovery of scattering information requires matrix inversions similar to DOT and spatial resolution may be ultimately limited by SNR and the requirement of using regularization. The present study may potentially lead to improved functional imaging, ability to visualize cancers and aid in tumor margin assessment (which often have enlarged cell nuclei and hence higher scattering coefficients). The improved quantitation could also lead to more confident mutliwavelength demixing for oxygen saturation estimation, and estimation of contrast agents or genetically encoded reporter distributions.