Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data

Peng Shao; Tyler Harrison; Roger J. Zemp

doi:10.1364/BOE.3.003240

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]

\frac{\partial Φ (r, t)}{\partial t} + c μ_{a} Φ (r, t) - c \nabla \cdot [D \nabla Φ (r, t)] = q (r, t)

where

Φ

is the optical fluence, 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 (μ_{a} + μ_{s}^{'})

. Here

μ_{a}

and

μ_{s}^{'}

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

μ_{s}^{'} ≫ μ_{a}

. With this hypothesis, fluence

Φ

will be almost isotropic. To calculate the fluence perturbations due to background heterogeneity at location

r

, we use the relations [16]

Φ (r) = Φ_{0} (r) + Φ_{S C} (r)

Φ_{S C} (r) = - \int \frac{δ μ_{a} (r')}{D_{0}} G_{0} (r, r') Φ (r') + \int \frac{δ D (r')}{D_{0}} \nabla G_{0} (r, r') \cdot \nabla Φ (r') d r'

Here

Φ_{0} (r)

and

Φ_{S C} (r)

are the background optical fluence and fluence perturbations, respectively.

G_{0} (r, r')

is the Green’s function representing propagation from

r \to r'

. Denote background diffusion coefficient and absorption coefficient by

D_{0}

and

μ_{a}

, the optical property distributions can be written as

D (r) = D_{0} + δ D (r)

and

μ_{a} (r) = μ_{a} + δ μ_{a} (r)

. The Born approximation is used in this paper to estimate fluence perturbations due to optical property variations at location

r

, when assuming that the fluence perturbations are small:

Φ_{S C} (r) ≪ Φ_{0} (r)

.

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

p^{o} = {[| p_{{111}}^{0} \dots p_{{T 11}}^{0} | \dots | p_{{1 M 1}}^{0} \dots p_{{T M 1}}^{0} | \dots | p_{{11 S}}^{0} \dots p_{{T 1 S}}^{0} | \dots | p_{{1 M S}}^{0} \dots p_{{T M S}}^{0} |]}^{T}

where

p_{{τ i k}}^{0}

is the observed pressure at time-point

τ

, with detector i, due to source k. The composite index

{τ i k}

can also be assigned the index

ρ

to map uniquely to the

ρ^{t h}

element of this vector such that

ρ = τ + (i - 1) T + (k - 1) M T

.

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{matrix} μ_{a} & D \end{matrix}]}^{T}$ , where $μ_{a}$ and D are rastorized column vectors of the optical absorption and diffusion coefficient maps. The objective is to find best estimate $\hat{u}$ for $u$ such that

\hat{u} = \underset{u}{\arg \min} ε (u)

where

ε (u) = {‖ p^{o} - p^{c} ‖}^{2} + λ β (u)

is the regularized least-squares error function with regularizing penalty function

β (u)

, which can be taken as a total variation function or other related function. 𝜆 is the so-called trade-off parameter that determines how much emphasis is placed on regularization while seeking the optimal solution. Forward model calculations of pressures received by transducer i due to source illumination k follow from the superposition of propagated initial pressure signals:

p_{{τ i k}}^{c} = \int \frac{p_{o, k} (r_{i}, t_{τ} - R_{i} / c)}{4 π R_{i}} d r = f^{k} (r_{i}, u, t_{τ})

with

R_{i} = | r_{i} - r |

where

r_{i}

is the location of the detector and

r

is a field point. This expression is valid for point-detectors with isotropic directivity and ideal electromechanical efficiency. Initial pressures are computed simply as

p_{o, k} (r, t) = Γ μ_{a} (r) Φ_{k} (r, u) δ (t),

where

Γ

is the Grüneisen parameter, assumed to be a constant in this paper, and

Φ_{k}

is the local laser fluence due to source k. To model the electromechanical response of the transducer,

f^{k} (r_{i}, u, t_{τ})

is convolved with impulse response

h_{r, x} (t)

. Thus in the expression for

p_{{τ i k}}^{c}

we can use a filtered version of the initial pressure distribution written as

p_{o, k} (r, t) = Γ μ_{a} (r) Φ_{k} (r, u) h_{r, x} (t) .

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 + δ u, t_{τ}) ≅ f^{k} (r_{i}, u, t_{τ}) + J_{i}^{k} δ u

, where

J

is the Jacobian matrix, and requiring that

\partial ε (u) / \partial δ u = 0

we arrive at the iterative algorithm

(J^{T} J) δ u = J^{T} (p^{o} - p^{c} (u)),

which can be written as

J^{T} J δ u = J^{T} Δ .

with

Δ \equiv p^{o} - p^{c} (u)

. For each iteration of the algorithm,

δ 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

δ 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

δ u = ξ H^{- 1} [J^{T} Δ]

, where

ξ

is a scalar chosen via a line-search algorithm to further minimize the error functional for each step. The Jacobian matrix has elements

J_{ρ j} = J_{{τ i k} j} = \frac{\partial f^{k} (r_{i}, u, t_{τ})}{\partial u_{j}} = \sum_{ℓ} (\frac{\partial f^{k} (r_{i}, u, t_{τ})}{\partial E_{ℓ}^{k}}) (\frac{\partial E_{ℓ}^{k}}{\partial u_{j}})

where

E_{ℓ}^{k} = μ_{a} (r_{ℓ}) Φ_{k} (r_{ℓ})

is the time-independent absorbed energy at location

r_{ℓ}

due to source k. It is straightforward to show that

{[α^{{τ i k}}]}_{ℓ} \equiv \frac{\partial f^{k} (r_{i}, u, t_{τ})}{\partial E_{ℓ}^{k}} = Γ \frac{h_{r x} (t_{τ} - R_{i ℓ} / c)}{4 π R_{i ℓ}}

which reduces to the Gruneisen parameter times the free-space acoustic Green’s function when

h_{r x} (t) \to δ (t)

and represents acoustic propagation from field point

ℓ

to detector i.

\partial E_{ℓ}^{k} / \partial u_{j}

represents optical absorbed energy variation with respect to optical property perturbations and it can be written as:

{[B^{k}]}_{j ℓ} \equiv \frac{\partial E_{ℓ}^{k}}{\partial μ_{j}} = {[\begin{matrix} \frac{\partial E_{ℓ}^{k}}{\partial μ_{a j}} & \frac{\partial E_{ℓ}^{k}}{\partial D_{j}} \end{matrix}]}^{T}

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

\frac{\partial E_{ℓ}^{k}}{\partial μ_{a j}} = δ_{ℓ j} Φ_{k} (r_{ℓ}) + μ_{a} (r_{ℓ}) \frac{\partial Φ_{k} (r_{ℓ})}{\partial μ_{a j}}

\frac{\partial E_{ℓ}^{k}}{\partial D_{j}} = μ_{a} (r_{ℓ}) \frac{\partial Φ_{k} (r_{ℓ})}{\partial D_{j}}

The partial derivatives can in turn be computed from the diffusion equation or in linearized form as

\frac{\partial Φ_{k} (r_{ℓ})}{\partial μ_{a j}} = - \frac{Δ x Δ y}{D_{0}} G_{0} (r_{ℓ}, r_{j}) Φ (r_{j}, r_{ℓ})

\frac{\partial Φ_{k} (r_{ℓ})}{\partial D_{j}} = - Δ x Δ y G_{0} (r_{ℓ}, r_{j}) \cdot Δ Φ (r_{j}, r_{ℓ})

Observed and computed pressure vectors have dimensions TMS × 1. For an N × N 2D grid of optical properties,

u

and

δ u

have dimensions 2N² × 1. Jacobian matrices have dimensions TMS × 2N². These matrices can be quite large. The Hessian H of size 2N² × 2N² has quadrant symmetry, reducing memory requirements. With the above notation, the

ρ t h

row of

J

is a 1 × 2N² vector given as

J_{P} = J_{{τ i k}} = {[B^{k} α^{{τ i k}}]}^{T}

where for a given k,

B^{k}

is a 2N² × 2N² matrix and for a given

{τ i k}

,

α^{{τ i k}}

is an 2N² × 1 column vector. The Jacobian needs not to be stored, but calculated row by row to find

b = J^{T} Δ = \sum_{ρ = 1}^{T M S} J_{ρ}^{T} J_{ρ}

and

H = \sum_{ρ = 1}^{T M S} J_{ρ}^{T} J_{ρ}

in a loop over

ρ

(a loop over the indices

{τ i k}

). Future Gradient-based optimization methods may circumvent the necessity of using such large matrices.

2.3. Inversion

Rewrite Eq. (7) as

H δ u = b

where

b = J^{T} Δ

. Now estimating the unknown parameters

δ u

is an inverse problem. Different techniques can be employed to solve the problem. For the i^th iteration, we minimize the following cost function

J_{i} (δ u_{i}) = {‖ H δ u_{i - 1} - b_{i - 1} ‖}_{2}^{2} + λ {‖ W δ u_{i - 1} ‖}_{2}^{2} .

For any vector

x

,

{‖ x ‖}_{2}^{2}

is the

ℓ_{2}

norm. The first term is the least squares norm which represents the misfit portion of the cost function. The second term is the Tickhonov regularization term and represents the constraints possessed on the parameter to be estimated. 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

δ {\hat{u}}_{i} = {(H_{i - 1}^{T} H_{i - 1} + λ W^{T} W)}^{- 1} H_{i}^{T} b .

Solution

δ {\hat{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 $Δ x = Δ y = 667 μ 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⁻¹ everywhere except a scattering perturbation of 1 cm⁻¹. The background absorption coefficient is taken as $μ_{a}$ = 0.1 cm⁻¹, 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_{r x} (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.

Fig. 1 The 2-dimensional simulation setup (a) and true optical property model (b), (c). Four light sources are distributed on each side of the imaging object, with a 3-mm gap backward from the soft tissue surface. 64 detectors are circumferentially located around the object.

Download Full Size | PDF

Fig. 2 (a) Normalized optical fluence distribution from source on top of the object. (b) Optical fluence perturbations due to only absorption heterogeneities. (c) Fluence perturbations due to only diffusion coefficient distribution. (d) Total fluence perturbations. (e) Ultrasound transducer channel data due to the same light source displayed in different forms.

Download Full Size | PDF

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 μ_a 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 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.

Fig. 3 Simulation results with the ratio metric method and the proposed method with 4 illuminations surrounding the object. From top to bottom, absorption and diffusion coefficient distribution obtained with the ratiometric method described in [15], and the 1st, 3rd, 5th, and 10th iteration with the proposed method.

Download Full Size | PDF

Fig. 4 Log of relative errors between true model and recovered absorption and diffusion coefficient distributions with different iterations.

Download Full Size | PDF

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.

Fig. 5 Singular values for different configurations. MIPAT with different source and detector count (1, 4 and 8 sources with 64 detectors and 16 detectors) are compared. Conditioning of CW-DOT with the same number of sources and detectors is also included for comparison.

Download Full Size | PDF

Table 1. Condition Number for Different Configurations

View Table

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 $ℓ_{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).

Case	Condition number
MI-PAT, 1 source, 64 detectors	3.7764e + 017
MI-PAT, 4 sources, 64 detectors	6.8934e + 007
MI-PAT, 4 sources, 16 detectors	1.3061e + 008
MI-PAT, 8 sources, 64 detectors	3.8724e + 007
CW-DOT, 8 sources	1.7838e + 013

Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data

Abstract

1. Introduction

2. Theory

2.1. Light propagation model

2.2. Reconstruction of the optical properties with ultrasound channel data

2.3. Inversion

3. Numerical simulation

4. Conclusions and discussion

Acknowledgments

References and links

Cited By

Figures (5)

Tables (1)

Equations (16)

Biomedical Optics Express