1. Introduction

Photoacoustic tomography holds great promise for quantitative optical imaging as it provides high quality images deeper in tissue than other optical contrast techniques. The images recovered directly from the data are actually estimates of the initial pressure distribution, which is a well-understood problem [1–3]. However quantitatively estimating optical absorption can be very difficult [4], as fluence is related to both optical absorption and scattering, and the Grüneisen parameter cannot be measured directly. In fact, it is been suggested by Bal and Ren that without multiple wavelengths, it is impossible to estimate more than two of the three contributing factors to initial pressure distribution (absorption, scattering, and the Grüneisen) uniquely [5].

Several attempts have been made at non-iterative solutions [6–8], but due to the interrelation between absorption and local fluence, iterative approaches have also shown some success [9, 10]. However, experimental work by Jetzfellner et al. [11] concluded that errors in optical scattering estimates can result in non-convergence. Multiple wavelengths have been used by Cox et al. [12] to recover absorption and scattering, and Bal and Ren [13] have proposed a method for also reconstructing the Grüneisen parameter.

Bal and Uhllmann posited that diffusion coefficients can be stably reconstructed from internal data (i.e. photoacoustic images) corresponding to a number of well-chosen boundary conditions (i.e. illuminations) [14]. Zemp introduced a ratiometric non-iterative method for reconstructing absorption distributions using multiple-illumination photoacoustic tomography (MIPAT) [15]. We have also previously discussed reconstruction of absorption, scattering, and Grüneisen distributions, and showed that while single illumination photoacoustic tomography is plagued by absorption-scattering non-uniqueness, multiple illumination photoacoustic tomography (MIPAT) can alleviate absorption scattering non-uniqueness and lead to a less ill-posed and better conditioned inverse problem [16]. Unfortunately, the ratiometric non-iterative approach used for quantitative MIPAT image reconstruction [15, 16] was susceptible to noise and applicable only for linearized problems in the diffusion regime with weak perturbations where the first Born approximation is valid. Additionally, iterative algorithms for reconstructing absorption and scattering distributions, suffer from computational complexity, numerical instability, and convergence problems connected with computation and inversion of very large and potentially ill-conditioned Jacobian and Hessian matrices, a problem that will become worse with larger datasets and with low signal-to-noise measurements. Shao et al. noted that non-ideal pressure reconstructions could lead to errors in quantitative photoacoustic tomography (PAT) and MIPAT estimates, then proposed an iterative algorithm that uses transducer channel data as a starting point [17]. Cox et al. have extended their previous multiwavelength work to include multiple illuminations as well [4]. Gao et al. introduced a Bregman method for multiple-illumination quantitative PAT [18]. Ren et al. have recently proposed a hybrid method using diffuse optical tomography to establish boundary values for the diffusion coefficient [19]. Still, for all approaches, the aforementioned problems of computational complexity and numerical stability are non-trivial, especially if surrounding optical properties are not well-known.

In this manuscript, we extend the fixed point iteration approach used by both Cox et al. [9] and Jetzfellner et al. [11] to multiple illuminations using an iterative least-squares approach. As was done in those works, we restrict our attention to reconstruction of absorption heterogeneities. Our proposed approach does not require inversion of large ill-conditioned Hessian matrices, is computationally simple, and robust to noise and inaccurate starting parameters, and offers substantially improved convergence properties compared to previous work. We investigate this technique at varying estimates of optical scattering and assess its stability in different noise conditions.

2. Theory

2.1. Single illumination

Assuming that an initial pressure distribution p₀ has already been reconstructed using any of the various techniques available [1–3], the problem for a single illumination is to reconstruct the optical absorption μ_a by modeling the fluence distribution Φ. The initial pressure distribution is actually a combination of μ_a, Φ and the Grüneisen parameter Γ which for simplicity we assume is uniform and constant. p₀ takes the form in equation 1.

Reconstruction is complicated by Φ, which should properly be written as Φ(r, μ_a, μ′_s): a function varying over spatial position r, absorption μ_a(r), and the reduced scattering coefficient μ′_s(r). For the simple iterative technique we are interested in [11], μ′_s is taken to be uniform. As in that work, we assume the diffusion equation is applicable, and use the finite element method solver available in the TOAST toolkit [20] to calculate Φ.

Assuming an initial guess of μ̂_a = 0, the iterative method of [9] proceeds as follows over iterations i: use the forward model to estimate Φⁱ; calculate an error parameter $\mathrm{\Delta }{p}_0=p_0-{\mu }_a^{(i)}{\mathrm{\Phi }}^{(i)}$; update ${\widehat{\mu }}_a^{(i+1)}=\frac{p_0}{{\mathrm{\Phi }}^{(i)}+\sigma }$ (where σ is a regularization parameter); and repeat with i = i + 1 until the error is within acceptable bounds, or the solution for μ̂_a has converged.

2.2. Extension to multiple illuminations

A similar derivation can be used for multiple illuminations. Consider that a number of detectors are used to reconstruct the initial pressure distributions ${\widehat{p}}_0^k\left(\mathit{r}\right)$ into an N × N image due to illuminations k = 1,...,S. The reconstructed initial pressures can be modeled as ${\widehat{p}}_0^k\left(\mathit{r}\right)=\mathrm{\Gamma }{\widehat{\mathrm{\Phi }}}_k\left(\mathit{r}\right){\widehat{\mu }}_a\left(\mathit{r}\right)$. Here, Φ̂_k(r) is the estimated fluence due to illumination k, and μ̂_a(r) is the estimated optical absorption coefficient distribution. For simplicity, Γ is considered to be spatially constant. From the M = N² pixels of each of the S reconstructed images, one can form a vector of the observation data ${\widehat{\mathit{p}}}_0-{\left[{\widehat{p}}_0^1\left({\mathit{r}}_1\right),\dots ,{\widehat{p}}_0^1\left({\mathit{r}}_M\right)\dots {\widehat{p}}_0^S\left({\mathit{r}}_1\right),\dots ,{\widehat{p}}_0^S\left({\mathit{r}}_M\right)\right]}^T$.

A vector equation may be formed to model computed initial pressure vectors: p̂₀ = Âμ̂_a, where μ̂_a = [μ̂_a(r₁),..., μ̂_a(r_M)]^T is an M × 1 column vector of estimated optical absorption coefficients, and where  = [Â₁,...,Â_S]^T with Â_k = Γ × diag(Φ̂_k), where Φ̂_k = [Φ̂_k(r₁),..., Φ̂_k(r_M)]^T.

The objective is to find μ̂_a such that the error between the observations and the computed images ε(μ̂_a) = ||p̂₀ − Âμ̂_a||² is a minimum. The least squares solution to this problem comes from solving (Â^TÂ)μ̂_a =Â^Tp̂₀ for the vector μ̂_a. Here ${\widehat{\mathit{A}}}^T\widehat{\mathit{A}}={\sum }_k{\widehat{\mathit{A}}}_k^T{\widehat{\mathit{A}}}_{\mathit{k}}={\mathrm{\Gamma }}^2\times \text{diag}\left({\sum }_k{\widehat{\mathrm{\Phi }}}_k^2\left({\mathit{r}}_1\right),\dots ,{\sum }_k{\widehat{\mathrm{\Phi }}}_k^2\left({\mathit{r}}_M\right)\right)$, and ${\widehat{\mathit{A}}}^T{\widehat{\mathit{p}}}_0={\sum }_k{\widehat{\mathit{A}}}_k^T{\widehat{\mathit{p}}}_0^k$ with ${\widehat{\mathit{p}}}_0=[{\widehat{p}}_0^k\left({\mathit{r}}_1\right),\dots ,{\widehat{p}}_0^k\left({\mathit{r}}_M\right)$. So, given the following form of the least-squares estimate: μ̂_a = (Â^TÂ)⁻¹Â^Tp̂₀, we have ${\widehat{\mathit{\mu }}}_{\mathit{a}}=\frac{1}{\mathrm{\Gamma }}\text{diag}\left(\frac{1}{{\sum }_k{\widehat{\mathrm{\Phi }}}_k^2\left({\mathit{r}}_1\right)},\dots ,\frac{1}{{\sum }_k{\widehat{\mathrm{\Phi }}}_k^2\left({\mathit{r}}_M\right)}\right)\mathit{B}$ with $\mathit{B}={\left[{\sum }_k{\widehat{\mathrm{\Phi }}}_k\left({\mathit{r}}_1\right){\widehat{p}}_0^k\left({\mathit{r}}_1\right),\dots ,{\sum }_k{\widehat{\mathrm{\Phi }}}_k\left({\mathit{r}}_M\right){\widehat{p}}_0^k\left({\mathit{r}}_M\right)\right]}^T$ and thus we obtain Eq. (2).

We can incorporate this into an iterative algorithm as follows: from initial fluence estimates Φ̂_k⁽⁰⁾(r, ${\widehat{\mathit{\mu }}}_a^{(0)}$) computed using a zeroth iteration approximation of optical absorption distribution ${\widehat{\mu }}_a^{(0)}$, we can compute a new estimate of optical absorption coefficients ${\widehat{\mu }}_a^{(1)}$ using the least-squares estimate above. This new estimate can be used to form a new fluence estimate using the diffusion equation or radiative transfer equation, and the process can be iterated until the error is sufficiently small. In order to avoid numerical instability, we introduce a regularization parameter β. β is intended only to ensure long-term convergence of the algorithm which may otherwise diverge due to noise effects, or in the case of experimental work, non-idealities. This results in Eq. (3) for iteration i + 1. For this equation, the absorption coefficients are guaranteed to be non-negative given that reconstructed initial pressures are themselves non-negative.

3. Simulations

Simulations were run in MATLAB using the TOAST [20] finite element-based forward solver for light propagation. In order to verify the simulation code, a similar absorption profile to that used in the Jetzfellner work [11] was simulated. In order to mimic the circular illumination used, 512 point sources located at one transport mean free path within the absorbing body were used. For the multiple illumination simulations, these sources were partitioned to give the appropriate number of images.

In order to measure the accuracy and convergence of MIPAT, we use a normalized root-mean-squared error (NRMSE) of the reconstructed absorption profile. NRMSE is calculated as: $NRMSE=\sqrt{\frac{\sum \sum {\left|{\widehat{\mu }}_a-{{\widehat{\mu }}_a}^{(i)}\right|}^2}{\sum \sum {\left|{\widehat{\mu }}_a\right|}^2}}$ where the summation is over the entire image, based on the ideally reconstructed absorption, ${\widehat{\mu }}_a=\frac{{\sum }_k{{\widehat{p}}_0}^k}{{\sum }_k{\mathrm{\Phi }}_k}$ (using the known Φ_k from the forward simulation), and μ̂_a⁽ⁱ⁾, the current quantitative image. This is the same measure as the second quality measure from the work of Jetzfellner et al. [11]. With this metric, we are interested in investigating three things: does increasing the number of sources improve convergence; what effect does the regularization parameter have; and how robust is the technique to noise. Convergence is not guaranteed to be a global minimum, only that the iteration will reach a stopping condition (i.e. further iteration no longer significantly changes the resulting reconstruction).

Our simulation of the experimental work provided very similar results using the same parameters used by Jetzfellner et al. [11] and a similar shaped phantom, shown in Fig. 1(a). That is, μ′_s = 20cm⁻¹, μ_a = 1.5cm⁻¹ for the inclusion, and μ_a = 0.2cm⁻¹ for the main body of the phantom. With σ = 0.001, we get very similar results as that work in terms of convergence with the varying estimates of μ′_s used for reconstruction. The most concerning aspect of the previous results is that the reconstruction does not converge even with the correct, uniform μ′_s = 20cm⁻¹. Figure 1(b) shows the initial pressure estimate and true fluence for a single uniform illumination. The first iteration of the algorithm gives a reasonable but inaccurate result seen in Fig. 1(c), but after 30 iterations, the solution is clearly diverging as in Fig. 1(d). We can instead use the four images with resulting fluences given in 1(e). By applying the MIPAT technique detailed above, setting β² = σ, and keeping all other parameters the same, we obtain very similar results for the first iteration in Fig. 1(f) (this is expected since it is the case with an initial guess of a uniform μ_a), and a much improved estimate of μ_a after 30 iterations in Fig. 1(g).

 

Fig. 1 (a) μ_a phantom; Single uniform illumination: (b) PA image (top) and true Φ (bottom), (c) first iteration, and (d) 30th iteration. Four illumination MIPAT: (e) PA images (top) and true Φ (bottom), (f) first iteration, and (g) 30th iteration. μ_a in cm⁻¹.

Download Full Size | PDF

Figure 2 shows more detailed results of the MIPAT simulation, comparing normalized root-mean-squared error (NRMSE) to the number of iterations for different numbers of illuminations. From this figure, we see that the previously non-convergent cases can be made to converge by increasing the number of illuminations using MIPAT techniques. Additionally, the speed of convergence improves with additional illuminations if not the absolute error.

 

Fig. 2 Simulated results of MIPAT with (a) 4, (b) 16, and (c) 512 illuminations with β = 0.032.

Download Full Size | PDF

Finally, we investigated the effect of noise on MIPAT reconstruction. As it turns out, β plays an important role in ensuring convergence in this case, so we explored the problem space in two dimensions: both over signal-to-noise ratio (SNR) and β. To make a fair comparison between different MIPAT illumination patterns, for each of S illuminations the total optical power delivered was $\frac{1}{S}$ the power used in the uniform illumination case (Fig. 1(b)) while noise levels per image remained the same. More precisely, the standard deviation of the noise added is based on the uniform illumination case at a level to produce the SNR defined as twenty times the base-10 logarithm of a ratio of signal (maximum signal in this case) to the standard deviation of the noise. This results in an SNR that is lower for each individual image in the MIPAT technique than the uniform illumination image. This should be somewhat equivalent to an experimental setup where portions of the delivered light are blocked to achieve each unique illumination. In cases with convergent solutions, 30 iterations was typically more than enough to demonstrate that behavior, so the root-mean-square error was measured at the 30th iteration for varying values of β and σ. Figure 3 shows our results with 4 and 16 illuminations.

 

Fig. 3 Simulated results of MIPAT after 30 iterations over different image signal-to-noise-ratios (SNR) and β values with (a) 4 and (b) 16 illuminations. The high SNR portion of (b) is presented in (c). Note that here, the standard deviation of the noise added to individual images is equal to that of the uniform illumination case.

Download Full Size | PDF

4. Discussion and conclusions

Figure 1(a–b) illustrates the problem that has been previously observed in experimental circumstances: even with a good estimate of μ′_s, a simple iterative technique does not necessarily result in a convergent solution. Our simulated results are in very good agreement with experiments from Jetzfellner et al. [11]. Non-convergence with good knowledge of background optical properties is highly problematic, and combined with the non-uniqueness problem that has previously been explored may prevent this technique from being practical.

MIPAT has already been demonstrated as a practical solution to the non-uniqueness problem [16], and our simulated results in Fig. 1 and Fig. 2 show how introducing multiple images from different optical excitations improves convergence - even where μ′_s is not well-known. While the speed of convergence is similar for the cases presented here, this figure shows that in all cases, the root-mean-squared error is minimized with the correct estimate for μ′_s as one would expect. Figure 2 illustrates an interesting behavior of the numerical method: for all values of μ′_s used in the reconstruction, a local minimum appears to be reached after two or three iterations before the error measure increases to eventually either diverge or converge to a different value. Indeed, fixed-point iteration numerical methods are not guaranteed to converge at a global minimum, and in fact are not guaranteed to converge at all.

While this improved convergence is an important result, the simulations in Fig. 1 and Fig. 2 do not include any noise or manipulation of the β regularization parameter. Figure 3 demonstrates the effects of SNR and β for 4 and 16 illuminations. In general, the trend that can be seen from these simulations is that increasing β can improve convergence in noisy situations, but that comes at the expense of reconstruction quality. The results seen in the figure are fairly intuitive, showing that low SNR hurts reconstruction quality, as does increasing β. However, there seems to be an ideal β for this set of simulations around 0.01 that provides reasonably accurate convergent solutions in low SNR conditions, while not significantly impacting reconstruction in high SNR conditions. While Fig. 3 demonstrates a range of SNR values, low SNR will be more typical of practical systems. In this case, more illuminations and larger regularization parameters are shown to be helpful. The choice of the regularization level using experimental data should be a topic of future work. A value of β that is too high will tend to cause the algorithm to converge in as little as a single iteration while providing inaccurate reconstruction. The choice of these regularization parameters (β for MIPAT and σ for the single illumination case) is certainly non-trivial, and in this work we were guided by the experimental work by Jetzfellner et al. [11], and we do find that the value of β that produced stable, yet accurate results is on the same order of magnitude as the σ = 0.001 used in that work if we use σ = β².

In comparing MIPAT with uniform illumination PAT, we chose to effectively divide the incident optical power into S sources each having $\frac{1}{S}$ of the available laser power. Even better performance might be expected if all energy of the incident laser were directed to each illumination location, however ANSI standards might limit fluence at each location in practice.

Our proposed algorithm makes no assumptions of simplified models of light transport or weak optical property perturbations. However, for simplicity, the light-propagation simulations used in this manuscript presently use the diffusion approximation. Like the work of Cox et al. [21], future work should assess robustness to cases where optical absorption coefficients are large or comparable to scattering coefficients and cases close to the point of entry, where traditional diffusion-regime approaches fail.

In the present form of the algorithm, reconstructed initial pressure distributions are taken as relatively faithful reconstructions of true initial pressures and it is assumed that there are known calibration factors relating reconstructed signals to true initial pressures. Additionally, no attempt was presently made to account for transducer spatio-temporal impulse-response or to otherwise account for non-ideal initial pressure reconstruction, and is a topic that should be given careful consideration in future work. Despite these limitations, faithful reconstructions were obtained on data which was extremely similar to the experimental data of Jetzfellner et al. [11], offering significant promise.

In this work, we have shown that MIPAT reconstruction can improve convergence over single illumination techniques, even where μ′_s is not well known. Not only that, but our investigation has shown that MIPAT can be used even in realistic noisy images by appropriately selecting the regularization parameter. While this technique does require a somewhat more complicated experimental setup than a single illumination method, the benefits of applying multiple illuminations may outweigh the drawbacks. Compared to previous multiple illumination algorithms, the fixed-point method discussed here does not require inversion of large Jacobian matrices [17,21] and hence is significantly more computationally efficient and stable.