Quantitative photoacoustic tomography augmented with surface light measurements

Olli Nykänen; Aki Pulkkinen; Tanja Tarvainen

doi:10.1364/BOE.8.004380

1. Introduction

Photoacoustic tomography (PAT) is an imaging modality based on the photoacoustic effect generated through the absorption of an externally introduced light pulse. The method combines optical contrast with high spatial resolution of ultrasound. The optical contrast is provided by distinctive absorption spectra by different chromophores. The chromophores of interest are, for example, haemoglobin, melanin and optionally various contrast agents. The ultrasonic waves carry this optical information directly to the surface with minimal scattering, thus retaining accurate spatial information as well. Nowadays, PAT can be used to provide images of soft biological tissues with high spatial resolution. It has successfully been applied to the visualisation of different structures in biological tissues, such as human blood vessels, microvasculature of tumors, and the cerebral cortex in small animals. For more information about PAT, see e.g. the reviews by [1–8] and the references therein.

Quantitative photoacoustic tomography (QPAT) is a technique which aims at estimating the absolute concentrations of the chromophores from photoacoustic images, i.e. from the reconstructed initial pressure distribution [9]. This is an ill-posed problem which needs to be approached in the framework of inverse problems. The concentrations of the chromophores can be estimated either by directly estimating them from photoacoustic images obtained at various wavelengths [10–16] or by first recovering the absorption coefficients at different wavelengths and then calculating the concentrations from the absorption spectra [9, 10, 13, 16]. Different approaches for the solution of the optical inverse problem of QPAT have been considered, see e.g. [17–35]. As an alternative to the conventional two stage approach, estimation of the optical parameters directly from the photoacoustic time-series has also been considered recently [34, 36–41]. In this work, the two stage approach is taken and it is assumed that the acoustic inverse problem, i.e. estimation of the initial pressure, has been solved. Further, one wavelength, i.e. estimation of absorption and scattering, is considered. Extensions of the developed numerical approach to one stage approach and spectral QPAT are straightforward.

It has been shown that, in order to obtain accurate estimates for absorption in QPAT, the scattering effects need to be taken into account [24,25,29,42]. Estimation of scattering is more ill-posed than absorption, and thus more sensitive to errors in data and modelling [25, 42]. Furthermore, it has been shown that false scattering values can lead to errors in absorption estimates [29].

In this work, improving QPAT by including information from surface measurements of light is investigated. Surface light measurements have been previously utilised in QPAT using the following approaches. In [43, 44], a two-step approach was suggested in which scattering distribution was first solved using diffuse optical tomography (DOT) measurements, and then this information was utilised in the estimation of the absorbed optical energy density in photoacoustic imaging. Furthermore, in [45], also utilising a two-step approach, a DOT experiment was first used to determine constant values for absorption and scattering, and then these were later utilised as background values in QPAT image reconstruction. In [46], a hybrid approach was introduced. In the approach, the vector field method developed in [24] was generalised and utilised in the estimation of the boundary parameters from surface light and the optical parameters inside the domain from absorbed optical energy.

In this work, estimating optical parameters using both absorbed optical energy density and surface light measurements is considered. In the approach, these two data sources are utilised simultaneously in order to solve the inverse problem of QPAT in a Bayesian framework. The approach is investigated in various imaging situations including differently sized domains and various internal structures of the optical parameters. Furthermore, in this paper, simultaneous estimation of the optical parameters and the Grüneisen coefficient is investigated. Simultaneous estimation of the absorption, scattering and Grüneisen parameters is non-unique if only one wavelength of light is used to obtain the QPAT data [13]. In this work, enhancing QPAT with surface light measurements is suggested to overcome this problem.

In QPAT, the most common approach has been to use the diffusion approximation (DA) as the light transport model, although the radiative transfer equation (RTE) has also been utilised [15, 21,25,38,47,48]. Recently, utilising Monte Carlo simulation methods have also been considered [35, 49]. In this work, the diffusion approximation (DA) is used as the forward model for light propagation. The DA is a valid approximation in a highly scattering medium further than a few scattering lengths from the light source, and thus it can be regarded as a relatively safe approximation for the simulation cases considered in this work.

The rest of the paper is organised as follows. The forward and inverse problem of QPAT are described in Section 2 and the related numerical implementations are given in Section 3. The results of simulations are shown in Section 4 and the conclusions are given in Section 5.

2. Methods

2.1. Forward model

In quantitative photoacoustic tomography, a short pulse of laser light is used to illuminate the tissue region of interest. The propagation of light can be described using the radiative transfer equation (RTE) [50]. In biomedical imaging, the RTE is often approximated with the diffusion approximation (DA)

\frac{1}{c} \frac{\partial Φ (r, t)}{\partial t} - \nabla \cdot \frac{1}{d (μ_{a} (r) + μ_{s}^{'} (r))} + \nabla Φ (r, t) + μ_{a} (r) Φ (r, t) = 0, r \in Ω

Φ (r, t) + \frac{1}{2 γ_{d}} \frac{1}{d (μ_{a} (r) + μ_{s}^{'} (r))} \frac{\partial Φ (r, t)}{\partial \hat{n}} = {\begin{array}{l} \frac{I_{s} (r, t)}{γ_{d}} & r \in ∊_{j} \\ 0, & r \in \partial Ω \ ∊_{j} \end{array}

where Φ(r, t) [W/mm²] is the photon fluence at time instance t, c is the speed of light in the medium, μ_a(r) [mm⁻¹] is the absorption coefficient, μ′_s(r) [mm⁻¹] is the (reduced) scattering coefficient, d is the dimension (d = 2, 3), I_s(r, t) [W/mm²] is a diffuse boundary current at the source position ∊_j ⊂ ∂Ω, ∂Ω is the boundary of the domain Ω, γ_d is a dimension-dependent constant which takes values γ₂ = 1/π and γ₃ = 1/4 and n̂ is an outward unit normal. The DA is a valid approximation in situations in which the radiance is almost an uniform distribution, i.e. in a scattering dominated medium further than a few scattering lengths from the light source [50].

As light propagates within the tissue, it is absorbed by chromophores. This generates localised increases in pressure. This pressure increase propagates through the tissue as an acoustic wave and is detected by ultrasound sensors at the surface of the tissue. The propagation of the acoustic wave occurs on a microsecond time scale, about five orders of magnitude slower than the optical propagation, so only the total absorbed optical energy density is of interest and not the rate of the absorption. Thus, in QPAT, light propagation can be modelled using a time-independent (continuous wave) model of light transport. The time-independent form of the DA is of the form

- \nabla \cdot \frac{1}{d (μ_{a} (r) + μ_{s}^{'} (r))} \nabla Φ (r) + μ_{a} (r) Φ (r) = 0, r \in Ω

Φ (r) + \frac{1}{2 γ_{d}} \frac{1}{d (μ_{a} (r) + μ_{s}^{'} (r))} \frac{\partial Φ (r)}{\partial \hat{n}} = {\begin{array}{l} \frac{I_{s} (r)}{γ_{d}}, & r \in ∊_{j} \\ 0, & r \in \partial Ω \ ∊_{j} \end{array}

where

Φ (r) = \int_{- \infty}^{\infty} Φ (r, t) d t

[J/mm²] is the time-independent fluence and I_s (r) [J/mm²] is a time-independent boundary light source. Furthermore, the absorbed optical energy density H(r) [J/mm³] is

H (r) = μ_{a} (r) Φ (μ_{a} (r), μ_{s}^{'} (r))

and the initial acoustic pressure p₀(r) [Pa] is [9]

p_{0} (r) = p (r, t = 0) = G (r) H (r)

where G(r) [unitless] is the Grüneisen parameter which is used to identify the photoacoustic efficiency. The time evolution of the resulting photoacoustic wave fields can be modelled using the equations of linear acoustics [51].

2.1.1. Measuring surface light

In this work, utilising surface light measurements in QPAT is investigated. Basically, this corresponds on using both QPAT and diffuse optical tomography data in determining the optical parameters. In a typical DOT measurement setup, the target is illuminated using either a short pulse of light (time-domain measurement systems), intensity modulated light (frequency-domain systems) or continuous light [52]. The measurable quantity is exitance Γ⁺(r,·) on the boundary of the target

Γ^{+} (r, \cdot) = - \frac{1}{d (μ_{a} (r) + μ_{s}^{'} (r))} \frac{\partial Φ (μ_{a} (r), μ_{s}^{'} (r), \cdot)}{\partial \hat{n}} = 2 γ_{n} Φ (μ_{a} (r), μ_{s}^{'} (r), \cdot) .

Depending on the measurement system, the exitance can be time-dependent Γ⁺(r, t) [W/mm²], frequency-dependent Γ⁺(r, ω) [J/mm²] or intensity only Γ⁺(r) [J/mm²]. The time-domain and frequency-domain data are related through Fourier-transform. In addition, other moments can be calculated from time-dependent measurements [53]. In this work, we consider frequency-domain data. In practice, these can be obtained either by a frequency-domain measurement device or by measuring the time-response of a time-domain system and Fourier-transforming the data into the frequency domain. Similarly, in order to solve the modelled exitance, the time-domain DA (1)–(2) or its counterpart in frequency domain

\frac{i ω}{c} Φ (r, ω) - \nabla \cdot \frac{1}{d (μ_{a} (r) + μ_{s}^{'} (r))} \nabla Φ (r, ω) + μ_{a} (r) Φ (r, ω) = 0, r \in Ω

Φ (r, ω) + \frac{1}{2 γ_{d}} \frac{1}{d (μ_{a} (r) + μ_{s}^{'} (r))} \frac{\partial Φ (r, ω)}{\partial \hat{n}} = {\begin{array}{l} \frac{I_{s} (r, ω)}{γ_{d}}, & r \in ∊_{j} \\ 0, & r \in \partial Ω \ ∊_{j} \end{array}

where ω is the angular modulation frequency of the input light and i is the imaginary unit, is used.

2.2. Inverse problem of QPAT

In this work, the optical inverse problem of QPAT, i.e. estimation of distributions of the optical parameters from photoacoustic images, is considered. The problem is approached in the framework of Bayesian inverse problems [28, 54]. A discrete observation model for QPAT in the presence of additive noise is

H_{meas} = H (μ_{a}, μ_{s}^{'}) + e_{qpat}

where H_meas ∈ ℝ^m_qpat is a measurement vector where m_qpat is the number of measurements which in this case is the number of illuminations multiplied with the number of discretised elements to represent the data space, μ_a = (μ_a₁, . . . , μ_{a_K})^T ∈ ℝ^K and μ′_s = (μ′_s₁, . . . , μ′_{s_K})^T ∈ ℝ^K are discretised absorption and scattering coefficients, K is the number of discretised parameters, H is the discretised forward operator corresponding to model (3)–(5) and e_qpat ∈ ℝ^m_qpat denotes the noise.

Let us assume that μ_a, μ′_s and H_meas are random variables. The solution of the inverse problem is the posterior probability density which is obtained through Bayesian formula and can be written in the form

π (μ_{a}, μ_{s}^{'} | H_{meas}) \propto π (H_{meas} | μ_{a}, μ_{s}^{'}) π (μ_{a}, μ_{s}^{'})

where π(H_meas|μ_a, μ′_s) is the likelihood density and π(μ_a, μ′_s) is the prior density. Assuming that the unknowns μ_a and μ′s and noise e are mutually independent and Gaussian distributed, i.e. μ_a ∼ 𝒩 (η_{μ_a}, Γ_{μ_a}),

μ_{s}^{'} ~ 𝒩 (η_{μ_{s}^{'}}, Γ_{μ_{s}^{'}})

and e ∼ 𝒩 (η_e,qpat, Γ_e,qpat), where η_{μ_a},

η_{μ_{s}^{'}}

and η_e,qpat are the means and Γ_{μ_a},

Γ_{μ_{s}^{'}}

and Γ_e,qpat are the covariances for the absorption, scattering and noise, respectively, the posterior density (11) becomes

\begin{array}{l} π (μ_{a}, μ_{s}^{'} | H_{meas}) \propto exp {- \frac{1}{2} {(H_{meas} - H (μ_{a}, μ_{s}^{'}) - η_{e, qpat})}^{T} Γ_{e, qpat}^{- 1} (H_{meas} - H (μ_{a}, μ_{s}^{'}) - η_{e, qpat}) \\ - \frac{1}{2} {(μ_{a}, η_{μ_{a}})}^{T} Γ_{μ_{a}}^{- 1} (μ_{a}, η_{μ_{a}}) - \frac{1}{2} {(μ_{s}^{'} - η_{μ_{s}^{'}})}^{T} Γ_{μ_{s}^{'}}^{- 1} (μ_{s}^{'} - η_{μ_{s}^{'}})} . \end{array}

The practical solution for the inverse problem is obtained by calculating point estimates from the posterior density. Since we are interested in computationally efficient inverse problem solvers, we consider here the maximum a posteriori (MAP) estimate. It is obtained as

\begin{array}{l} ({\hat{μ}}_{a}, {\hat{μ}}_{s}^{'}) & = \underset{(μ_{a}, μ_{s}^{'})}{arg max} π (μ_{a}, μ_{s}^{'} | H_{meas}) \\ = \underset{(μ_{a}, μ_{s}^{'})}{arg min} {\frac{1}{2} {‖ L_{e, qpat} (H_{meas} - H (μ_{a}, μ_{s}^{'})) - η_{e, qpat} ‖}_{2}^{2} \\ + \frac{1}{2} {‖ L_{μ_{a}} (μ_{a} - η_{μ_{a}}) ‖}_{2}^{2} + \frac{1}{2} {‖ L_{μ_{s}^{'}} (μ_{s}^{'} - η_{μ_{s}^{'}}) ‖}_{2}^{2}} \end{array}

where L_e,qpat is the Cholesky decomposition of the inverse of the noise covariance matrix

Γ_{e, qpat}^{- 1} = L_{e, qpat}^{T} L_{e, qpat}

, and L_{μ_a} and

L_{μ_{s}^{'}}

are the Cholesky decompositions of the inverse of the prior covariance matrices for absorption and scattering,

Γ_{μ_{a}}^{- 1} = L_{μ_{a}}^{T} L_{μ_{a}}

and

Γ_{μ_{s}^{'}}^{- 1} = L_{μ_{s}^{'}}^{T} L_{μ_{s}^{'}}

, respectively. Thus, in the image reconstruction problem of QPAT, we seek to find the discretised distributions of absorption and scattering coefficients (μ̂_a, μ̂′_s) which solve the minimisation problem (13).

2.2.1. QPAT augmented with surface light data

A discrete observation model for DOT in the presence of additive noise is

Γ_{meas}^{+} = Γ^{+} (μ_{a}, μ_{s}^{'}) + e_{dot}

where

Γ_{meas}^{+} \in ℝ^{m_{dot}}

is a measurement vector of light exitance measured at the detectors on the surface of the target where m_dot is the number of measurement, Γ⁺ is the discretised forward operator corresponding to model (7)–(9) and e_dot ∈ ℝ^m_dot denotes the noise. Following a similar derivation as for QPAT, the inverse problem of QPAT augmented with surface light measurements can formulated as a minimisation problem

\begin{array}{r} ({\hat{μ}}_{a}, {\hat{μ}}_{s}^{'}) = \underset{(μ_{a}, μ_{s}^{'})}{arg min} {\frac{1}{2} {‖ L_{e, qpat} (H_{meas} - H (μ_{a}, μ_{s}^{'}) - η_{e, qpat}) ‖}_{2}^{2} \\ + \frac{1}{2} {‖ L_{e, dot} (Γ_{meas}^{+} - Γ^{+} (μ_{a}, μ_{s}^{'}) - η_{e, dot}) ‖}_{2}^{2} \\ + \frac{1}{2} {‖ L_{μ_{a}} (μ_{a}, η_{μ_{a}}) ‖}_{2}^{2} + {‖ L_{μ_{s}^{'}} (μ_{s}^{'}, η_{μ_{s}^{'}}) ‖}_{2}^{2}} \end{array}

where η_e,dot is the mean of the noise of the surface light measurements and L_e,dot is the Cholesky decomposition of the inverse of the noise covariance matrix of the surface light measurements

Γ_{e, dot}^{- 1} = L_{e, dot}^{T} L_{e, dot}

.

2.2.2. Simultaneous estimation of the Grüneisen parameter

Generally, estimation of the absorption, scattering and Grüneisen parameter is non-unique if only one wavelength of light is used [24]. This has been overcome by using multiple optical wavelengths and estimating the optical parameters and the Grüneisen coefficient simultaneously [13, 15, 16]. In this work, simultaneous estimation of the absorption, scattering and Grüneisen parameter is considered using initial pressure and surface light as data. The minimisation problem is of the form

\begin{array}{l} ({\hat{μ}}_{a}, {\hat{μ}}_{s}^{'}, \hat{G}) = & \underset{(μ_{a}, μ_{s}^{'}, G)}{arg min} {\frac{1}{2} {‖ L_{e, p} (p_{0, meas} - p_{0} (μ_{a}, μ_{s}^{'}, G)) - η_{e, p} ‖}_{2}^{2} \\ + \frac{1}{2} {‖ L_{e, dot} (Γ_{meas}^{+} - Γ^{+} (μ_{a}, μ_{s}^{'}) - η_{e, dot}) ‖}_{2}^{2} + \frac{1}{2} {‖ L_{μ_{a}} (μ_{a} - η_{μ_{a}}) ‖}_{2}^{2} \\ + \frac{1}{2} {‖ L_{μ_{s}^{'}} (μ_{s}^{'}, η_{μ_{s}^{'}}) ‖}_{2}^{2} + \frac{1}{2} {‖ L_{G} (G - η_{G}) ‖}_{2}^{2}} \end{array}

where G = (G₁, . . . , G_K)^T ∈ ℝ^K is the discretised distribution of the Grüneisen parameter which was assumed to be Gaussian distributed and independent of the noise, p_0,meas is the initial pressure obtained from measurements, p₀(μ_a, μ′_s, G) is the modelled initial pressure, and η_e,p is the mean and L_e,p is the Cholesky decomposition of the inverse of the noise covariance matrix

Γ_{e, p}^{- 1} = L_{e, p}^{T} L_{e, p}

of the initial pressure data. Further, η_G is the mean and L_G is the Cholesky decomposition of the inverse of the covariance matrix

Γ_{G}^{- 1} = L_{G}^{T} L_{G}

of the prior for the Grüneisen parameter.

3. Numerical implementations

3.1. FE-approximation of the DA

In this work, a finite element method (FEM) is used for the numerical solution of the DA. In order to obtain the FE-approximation, first a variational problem is formulated of the DA with its boundary condition, and then the solution of the variational problem is approximated in a piece-wise linear basis. As a result, a matrix equation can be derived. For more details, see e.g. [25,55].

In the case of QPAT, the FE-approximation of the time-independent DA with the diffuse boundary source model, Eqs. (3)–(4), is obtained by solving the matrix equation

A_{qpat} Φ_{qpat} = b_{qpat}

where Φ_qpat is the fluence in the nodes of the FE-mesh and A_qpat = K + C + R. In the case of the frequency-domain DA, Eqs. (8)–(9), the FE-approximation of the DA on an angular modulation frequency ω is obtained by solving the matrix equation

A_{dot} Φ_{dot} = b_{dot}

where Φ_dot is the fluence in the nodes of the FE-mesh and A_dot = K + C + R + Z. The matrices K, C, R and Z and the source vectors b_qpat and b_dot are of the form

K (p, q) = \int_{Ω} \frac{1}{d (μ_{a}, μ_{s}^{'})} \nabla φ_{q} (r) \cdot \nabla φ_{p} (r) d r

C (p, q) = \int_{Ω} μ_{a} φ_{q} (r) φ_{p} (r) d r

R (p, q) = \int_{\partial Ω} 2 γ_{d} φ_{q} (r) φ_{p} (r) d S

Z (p, q) = \frac{i ω}{c} \int_{Ω} φ_{q} (r) φ_{p} (r) d r

b_{qpat} (p) = \int_{ε_{j}} 2 I_{s} φ_{p} (r) d S

b_{dot} (p) = \int_{ε_{j}} 2 I_{s} (ω) φ_{p} (r) d S

where φ_q (r) and φ_p (r) are FE-basis functions, q, p = 1, . . . , N, and N is the number of nodes.

The QPAT data, absorbed optical energy density, H can be computed from the fluence Φ_qpat obtained with (17) using Eq. (5) and the surface light measurements Γ⁺ can be computed from the fluence Φ_dot obtained with (18) using Eq. (7).

3.2. Gauss-Newton iteration

In this work, the minimisation problems (13), (15) and (16) are solved using a Gauss-Newton method. In the case of QPAT image reconstruction problem augmented with surface light measurements, Eq. (15), the Gauss-Newton iteration can be written in a form

x_{(i + 1)} = x_{(i)} + s_{(i)} {(J_{(i)}^{T} L_{e}^{T} L_{e} J_{(i)} + L_{x}^{T} L_{x})}^{- 1} (J_{(i)}^{T} L_{e}^{T} L_{e} (F_{meas} - F_{(i)} - η_{e}) - L_{x}^{T} L_{x} (x_{(i)} - η_{x}))

where x = (μ_a, μ′_s)^T = (μ_a₁ , . . . , μ_{a_K}, μ′_s₁, . . . , μ′_{s_K})^T ∈ ℝ^2K are the estimated absorption and scattering parameters which in this work are represented in piece-wise constant bases

μ_{a} (r) \approx \sum_{k = 1}^{K} μ_{a_{k}} χ_{k}^{(μ_{a})} (r)

and

μ_{s}^{'} (r) \approx \sum_{k = 1}^{K} μ_{s_{k}}^{'} χ_{k}^{(μ_{s}^{'})} (r)

where χ_k (r) is a characteristic function of the element Ω_k. It should be noted that other bases can also be used for the representation of the optical parameters and that different discretisations can be used for different parameters. Furthermore,

F_{meas} = {(H_{meas}, Γ_{meas}^{+})}^{T}

are the measured absorbed optical energy density and exitance, F =(H, Γ⁺)^T is the solution of the discretised forward models for QPAT, Eqs. (17) and (5), and DOT, Eqs. (18) and (7), s is the step length of the minimisation algorithm determined using a line-search method, and

L_{e} = (\begin{matrix} L_{e, qpat} & 0 \\ 0 & L_{e, dot} \end{matrix}), η_{e} = (\begin{array}{l} η_{e, qpat} \\ η_{e, dot} \end{array}), L_{x} = (\begin{matrix} L_{μ_{a}} & 0 \\ 0 & L_{μ_{s}^{'}} \end{matrix}), η_{x} = (\begin{array}{l} η_{μ_{a}} \\ η_{μ_{s}^{'}} \end{array}) .

Further, J is the Jacobian of the form

J = (\begin{matrix} J_{μ_{a}}^{qpat} & J_{μ_{s}^{'}}^{qpat} \\ J_{μ_{a}}^{dot} & J_{μ_{s}^{'}}^{dot} \end{matrix})

where Jacobians for absorption and scattering are

J_{μ_{a}}^{qpat} = (j_{μ_{a}}^{{qpat}^{(1)}}, \dots, j_{μ_{a}}^{{qpat}^{(K)}})

,

J_{μ_{s}^{'}}^{qpat} = (j_{μ_{s}^{'}}^{{qpat}^{(1)}}, \dots, j_{μ_{s}^{'}}^{{qpat}^{(K)}})

,

J_{μ_{a}}^{dot} = (j_{μ_{a}}^{{dot}^{(1)}}, \dots, j_{μ_{a}}^{{dot}^{(K)}})

and

J_{μ_{s}^{'}}^{dot} = (j_{μ_{s}^{'}}^{{dot}^{(1)}}, \dots, j_{μ_{s}^{'}}^{{dot}^{(K)}})

where the columns k = 1, . . . , K of each matrix are obtained by vectorisation of

j_{μ_{a}}^{{qpat}^{(k)}} = - μ_{a}^{(k)} ℳ_{qpat} A_{qpat}^{- 1} \frac{\partial A_{qpat}}{\partial μ_{a_{k}}} A_{qpat}^{- 1} b_{qpat} + χ_{k} ℳ_{qpat} A_{qpat}^{- 1} b_{qpat}

j_{μ_{a}}^{{dot}^{(k)}} = - ℳ_{dot} A_{dot}^{- 1} \frac{\partial A_{dot}}{\partial μ_{a_{k}}} A_{dot}^{- 1} b_{dot}

j_{μ_{s}^{'}}^{{qpat}^{(k)}} = - μ_{a}^{(k)} ℳ_{qpat} A_{qpat}^{- 1} \frac{\partial A_{qpat}}{\partial μ_{s_{k}}^{'}} A_{qpat}^{- 1} b_{qpat}

j_{μ_{s}^{'}}^{{dot}^{(k)}} = - ℳ_{dot} A_{dot}^{- 1} \frac{\partial A_{dot}}{\partial μ_{s_{k}}^{'}} A_{dot}^{- 1} b_{dot}

where

\frac{\partial A_{qpat}}{\partial μ_{a_{k}}} (p, q) = \frac{\partial A_{dot}}{\partial μ_{a_{k}}} (p, q) = - \frac{1}{d {(μ_{a_{k}} + μ_{s_{k}}^{'})}^{2}} \int_{Ω_{k}} \nabla φ_{q} (r) \cdot \nabla φ_{p} (r) d r + \int_{Ω_{k}} φ_{q} (r) φ_{p} (r) d r

\frac{\partial A_{qpat}}{\partial μ_{s_{k}}^{'}} (p, q) = \frac{\partial A_{dot}}{\partial μ_{s_{k}}^{'}} (p, q) = - \frac{1}{d {(μ_{a_{k}} + μ_{s_{k}}^{'})}^{2}} \int_{Ω_{k}} \nabla φ_{q} (r) \cdot \nabla φ_{p} (r) d r

and ℳ_qpat ∈ ℝ^m_qpat×N is a measurement matrix which contains the discretised measurement operator for QPAT which maps the piece-wise linear solution of the QPAT forward model to piece-wise constant data space and ℳ_dot ∈ ℝ^m_dot×N is a measurement matrix which contains the discretised measurement operator for surface light which maps the piece-wise linear solution of the DOT forward model to detector measurements.

The Gauss-Newton algorithm for the solution of the conventional QPAT problem (13), can be formulated from (25) by dropping the DOT related matrices and vectors from the forward model F, measurement vector F_meas, noise statistics L_e, η_e and the Jacobian matrix J.

In the case of estimating the absorption, scattering and Grüneisen parameter simultaneously from the initial pressure and surface light data, i.e. the minimisation problem (16), the Gauss-Newton iteration is of the form

{\tilde{x}}_{(i + 1)} = {\tilde{x}}_{(i)} + s_{(i)} {({\tilde{J}}_{(i)}^{T} {\tilde{L}}_{e}^{T} {\tilde{L}}_{T} {\tilde{J}}_{(i)} + {\tilde{L}}_{x}^{T} {\tilde{L}}_{x})}^{- 1} ({\tilde{J}}_{(i)}^{T} {\tilde{L}}_{e}^{T} {\tilde{L}}_{e} ({\tilde{F}}_{meas} - {\tilde{F}}_{(i)} - {\tilde{η}}_{e}) - {\tilde{L}}_{x}^{T} {\tilde{L}}_{x} ({\tilde{x}}_{(i)} - η_{\tilde{x}}))

where where x̃ = (μ_a, μ′_s, G)^T = (μ_a₁ , . . . , μ_{a_K}, μ′_s₁, . . . , μ′_{s_K}, G₁, . . . , G_K)^T ∈ ℝ^3K are the estimated absorption, scattering and Grüneisen parameters,

{\tilde{F}}_{meas} = {(p_{0, meas}, Γ_{meas}^{+})}^{T}

are the measured initial pressure and exitance, F̃_meas = (p₀, Γ⁺)^T is the solution of the discretised forward models and

{\tilde{L}}_{e} = (\begin{matrix} L_{e, p} & 0 \\ 0 & L_{e, dot} \end{matrix}), {\tilde{η}}_{e} = (\begin{matrix} η_{e, p} \\ η_{e, dot} \end{matrix}), {\tilde{L}}_{x} = (\begin{matrix} L_{μ_{a}} & 0 & 0 \\ 0 & L_{μ_{s}^{'}} & 0 \\ 0 & 0 & L_{G} \end{matrix}), η_{\tilde{x}} = (\begin{matrix} η_{μ_{a}} \\ η_{μ_{s}^{'}} \\ η_{G} \end{matrix}) .

Further, J̃ is the Jacobian of the form

\tilde{J} = (\begin{matrix} {\tilde{J}}_{μ_{a}}^{qpat} & {\tilde{J}}_{μ_{s}^{'}}^{qpat} & {\tilde{J}}_{G}^{qpat} \\ J_{μ_{a}}^{dot} & J_{μ_{s}^{'}}^{dot} & 0 \end{matrix})

where Jacobians for the absorption, scattering and Grüneisen parameter are

{\tilde{J}}_{μ_{a}}^{qpat} = ({\tilde{j}}_{μ_{a}}^{{qpat}^{(1)}}, \dots, {\tilde{j}}_{μ_{a}}^{{qpat}^{(K)}})

,

{\tilde{J}}_{μ_{s}^{'}}^{qpat} = ({\tilde{j}}_{μ_{s}^{'}}^{{qpat}^{(1)}}, \dots, {\tilde{j}}_{μ_{s}^{'}}^{{qpat}^{(K)}})

and

{\tilde{J}}_{G}^{jqpat} = ({\tilde{j}}_{G}^{{qpat}^{(1)}}, \dots, {\tilde{j}}_{G}^{{qpat}^{(K)}})

where the columns k = 1, . . . , K of each matrix are obtained by vectorisation of

{\tilde{j}}_{μ_{a}}^{{qpat}^{(k)}} = - G^{(k)} μ_{a}^{(k)} ℳ_{qpat} A_{qpat}^{- 1} \frac{\partial A_{qpat}}{\partial μ_{a_{k}}} A_{qpat}^{- 1} b_{qpat} + G^{(k)} ℳ_{qpat} A_{qpat}^{- 1} b_{qpat}

{\tilde{j}}_{μ_{s}^{'}}^{{qpat}^{(k)}} = - G^{(k)} μ_{a}^{(k)} ℳ_{qpat} A_{qpat}^{- 1} \frac{\partial A_{qpat}}{\partial μ_{s_{k}}^{'}} A_{qpat}^{- 1} b_{qpat}

{\tilde{j}}_{G}^{{qpat}^{(k)}} = - μ_{a}^{(k)} ℳ_{qpat} A_{qpat}^{- 1} b_{qpat}

and

J_{μ_{a}}^{dot}

and

J_{μ_{s}^{'}}^{dot}

are as in (28) and (30).

4. Results

The proposed approach was tested with simulations which were implemented in two-dimensional rectangular domains. Estimating absorption and scattering was investigated with varying domain size and varying parameter distributions using absorbed optical energy density and exitance as data. Furthermore, simultaneous estimation of the absorption, scattering and Grüneisen parameter was investigated using initial pressure and exitance as data. The results were compared to conventional QPAT reconstructions. The simulations were performed with a Windows workstation with Intel(R) processor with 2 cores, a speed of 3.16 GHz and 8 GB RAM using matlab (R2014a, The MathWorks Inc. Natick, Massachusetts, USA).

4.1. Data simulation

In all of the simulations, four planar illuminations, one at each side of the rectangle, with a uniform intensity covering the whole side length and a total energy of 1 J, were used. It should be noted that, since the numerical simulations are performed with a noise amplitude related to the data, the absolute magnitude of the input light energy does not affect the results of these simulations. The absorbed optical energy density, i.e. QPAT data, was simulated by using the FE-solution of the DA (17) to obtain the fluence and then multiplying that with the absorption to get the absorbed optical energy density using Eq. (5). To get the initial pressure distribution, the absorbed optical energy density was multiplied with the Grüneisen parameter using Eq. (6). In all simulations, Gaussian distributed noise with zero mean and standard deviation of 1 % of the maximum absorbed optical energy density or initial pressure was added to the simulated data.

Furthermore, the boundary light data was simulated by first solving the FE-approximation of the frequency domain DA (18) for fluence, and then calculating the exitance using Eq. (7). The exitance was solved in 174 equally distributed detector points (58 at each side) on the object boundary excluding the illumination side of the target. In the simulations, the angular modulation frequency of ω = 100·10⁶ rad/s was used. Gaussian distributed noise with standard deviation of 1 % of the amplitude and phase was added to the simulated data.

The number of nodes and elements of the FE-discretisations used in the simulations are given in Table 1. The optical parameters were represented in piece-wise constant bases using the elements of the FE-discretisation.

Table 1. The number of nodes and elements in the FE-discretisations used in simulating the data.

View Table | View all tables in this article

4.2. Reconstructing the parameters of interest

The absorption and scattering distributions were reconstructed using the suggested augmented QPAT approach by minimising (15). The minimisation problem was solved using the Gauss-Newton method (25) equipped with a line search algorithm for the determination of the step length. The results were compared to the conventional QPAT reconstructions by minimising (13). Furthermore, the absorption, scattering and Grüneisen parameter distributions were reconstructed using the suggested augmented QPAT approach by minimising (16). The minimisation problem was solved using the Gauss-Newton method (33) equipped with a line search algorithm for the determination of the step length. All minimisation problems converged in less than 15 iterations. Typically, the augmented QPAT took few more simulations to converge than the conventional QPAT. The number of nodes and elements of the FE-discretisations used in the reconstructions is given in Table 2. The optical parameters were represented in piece-wise constant bases by the elements of the FE-discretisations.

Table 2. The number of nodes and elements in the FE-discretisations used in the reconstructions.

View Table | View all tables in this article

In this work, the prior model for the unknown parameters μ_a, μ′_s and Γ was chosen to be based on the Ornstein-Uhlenbeck process [16, 56]. The prior is a Gaussian distribution with mean η and covariance matrix Γ which was defined as being proportional to

Γ = σ^{2} Ξ

where σ is the standard deviation of the prior and Ξ is a matrix which has its elements defined as

Ξ_{i j} = exp (- ‖ r_{i} - r_{j} ‖ / ξ),

where i and j denote the row and column indices of the matrix, r_i and r_j denote the coordinates of the centre of elements i and j, and ξ is the characteristic length scale of the prior describing the spatial distance that the parameter is expected to have (significant) spatial correlation for. In this work, the mean and standard deviation of the prior were chosen based on expected mean and variation that the parameters were assumed to have, and the characteristic length scale was chosen based on size of the structures that the medium was assumed to have. The mean, standard deviation and characteristic length scale used in the simulations are given in Table 3.

Table 3. The mean and standard deviation of the prior distributions for the absorption η_{μ_a} (mm⁻¹), σ_{μ_a} (mm⁻¹), scattering $η_{μ_{s}^{'}}$ (mm⁻¹), $σ_{μ_{s}^{'}}$ (mm⁻¹) and Grüneisen parameter η_G (unitless), σ_G (unitless), and the characteristic length scale ξ(mm) used in the simulations.

View Table | View all tables in this article

The validity of the reconstructions was evaluated by computing the mean relative errors for the estimates

E_{μ_{a}} = 100 % \cdot \frac{{‖ {\hat{μ}}_{a} - μ_{a}^{TRUE} ‖}_{2}}{{‖ μ_{a}^{TRUE} ‖}_{2}}

E_{μ_{s}^{'}} = 100 % \cdot \frac{{‖ μ_{s}^{'} - μ_{s}^{' TRUE} ‖}_{2}}{{‖ μ_{s}^{' TRUE} ‖}_{2}}

E_{G} = 100 % \cdot \frac{{‖ \hat{G} - G^{TRUE} ‖}_{2}}{{‖ G^{TRUE} ‖}_{2}}

where

μ_{a}^{TRUE}

,

μ_{s}^{' TRUE}

and G^TRUE are the simulated parameters interpolated to the reconstruction grid and μ̂_a, μ̂′_s and Ĝ are the estimated parameters.

4.3. Results

4.3.1. Varying domain size

The performance of the augmented QPAT in different size domains was investigated. The sizes of the simulation domains were 20 mm × 20 mm, 40 mm × 40 mm and 60 mm × 60 mm. The absorption and scattering inclusions were stripe-like structures with thickness of 2 mm and length of the side-length of the domain minus 1 mm. The stripes were located with a distance of 3 mm from each other. The simulated parameter distributions and the reconstructions are shown in Fig. 1. The relative errors for the estimated absorption and scattering, Eqs. (40)–(41), are given in Table 4.

Fig. 1 Absorption (columns 1–3) and scattering (columns 4–6) distributions in different size domains: 20 mm ×20 mm (first row), 40 mm×40 mm (second row) and 60 mm×60 mm (third row). Columns from left to right: simulated absorption (first column), reconstructed absorption obtained using the augmented (second column) and conventional (third column) QPAT, simulated scattering (fourth column), reconstructed scattering obtained using augmented (fifth column) and conventional (sixth column) QPAT.

Download Full Size | PDF

Table 4. Mean relative errors for absorption E_{μ_a} (%) and scattering $E_{μ_{s}^{'}}$ (%) obtained with augmented QPAT and conventional QPAT.

View Table | View all tables in this article

As it can been seen, the absorption distributions reconstructed using the augmented QPAT and the conventional QPAT approaches are almost the same quality in the 20 mm × 20 mm size domain. The scattering estimates, on the other hand, obtained with the augmented QPAT resemble the original distribution more than the estimates obtained with the conventional QPAT. The difference between the two approaches becomes more apparent when the domain size increases. Especially in the domain of size 40 mm × 40 mm, the reconstructions are clearly better quality when the surface light data has been utilised in the reconstructions. In the largest domain, 60 mm × 60 mm, the quality of both conventional and augmented QPAT is weaker than in the smaller domains. That is, one is reaching the limit where QPAT reconstructions cannot significantly be improved by including exitance data.

The calculated mean relative errors, Table 4, support these results. The relative errors of the absorption and scattering estimates obtained with the augmented QPAT are smaller than those of the conventional QPAT reconstructions. The most accurate estimates are obtained in the 20 mm × 20 mm domain and as the domain size increases the relative errors of the estimates increase as well. The most significant improvement into the accuracy of the estimates by the augmented QPAT is obtained in 40 mm × 40 mm for absorption estimates.

The convergence of the minimisation problems was visualised by calculating the mean relative errors for absorption and scattering, Eqs. (40)–(41), at each iteration and plotting the results. All the simulations had a similar behaviour, and thus only the results of once case (simulations in 20 mm × 20 mm domain) are shown in Fig. 2. As it can be seen, the augment QPAT required more iterations to converge than the conventional QPAT but it converged to more accurate estimates for absorption and scattering.

Fig. 2 Mean relative errors for absorption E_{μ_a} (left image) and scattering $E_{μ_{s}^{'}}$ (right image) against iteration obtained with augmented QPAT and conventional QPAT in 20 mm ×20 mm domain.

Download Full Size | PDF

4.3.2. Variations in the optical parameter distribution

The capability of the augmented QPAT approach to reconstruct fine structures in internal optical parameter distributions was investigated. The simulation domain size was 20 mm ×20 mm. In all simulations, the absorption and scattering were given two different values but their distribution was varied from coarse to fine structures. The total area of absorbing and scattering inclusions was kept the same in all cases. The simulated and the reconstructed absorption and scattering distributions are shown in Fig. 3. The mean relative errors for absorption and scattering, Eqs. (40)–(41), are given in Table 4.

Fig. 3 Absorption (columns 1–3) and scattering (columns 4–6) parameters with different distributions of the optical parameters (rows 1–4). Columns from left to right: simulated absorption (first column), reconstructed absorption obtained using the augmented (second column) and conventional (third column) QPAT, simulated scattering (fourth column), reconstructed scattering obtained using the augmented (fifth column) and conventional (sixth column) QPAT.

Download Full Size | PDF

As it can be seen, QPAT augmented with surface light measurements gives better quality reconstructions both for absorption and scattering. In the absorption reconstructions, the shapes of the absorption inclusions are well defined by both methods. However, it seems that the augmented QPAT gives more accurate values for the absorbing inclusions. In the case of scattering, the inclusions are better distinguished if augmented QPAT data is used. The relative errors of the absorption and scattering estimates obtained with the augmented QPAT are smaller than those obtained with the conventional QPAT. That is, the augmented QPAT provides more accurate quantitative estimates for the optical parameters.

4.3.3. Simultaneous estimation of the Grüneisen parameter

The simultaneous estimation of the absorption, scattering and Grüneisen parameter was investigated in 10 mm × 10 mm domain. For comparison, two types of conventional QPAT reconstructions, Eq. (13), were computed. In the first one, the Grüneisen parameter was assumed to have a constant value G = 0.3, which is the mean of the simulated Grüneisen parameter distribution. In the second one, the simulated Grüneisen parameter distribution was mapped to the reconstruction mesh. This can be considered the best possible reference for the reconstructions, which however is unrealistic since, in practice, the Grüneisen parameter distribution is unknown. The simulated and reconstructed parameter distributions are shown in Fig. 4. The mean relative errors, Eqs. (40)–(42), are given in Table 5.

Fig. 4 Absorption (first row), scattering (second row) and Grüneisen (third row) parameters. Columns from left to right: simulated distributions (first column), reconstructions obtained using the augmented QPAT (second column), and reconstructions obtained using the conventional QPAT with a constant (third column) and correct (fourth column) values of the Grüneisen parameter.

Download Full Size | PDF

Table 5. Mean relative errors for absorption E_{μ_a} (%), scattering $E_{μ_{s}^{'}}$ (%) and Grüneisen parameter E_G (%) obtained with augmented QPAT and conventional QPAT.

View Table | View all tables in this article

As it can be seen, the approach can produce as good quality reconstructions for absorption and scattering as the reference approach with the correct Grüneisen parameter distribution and better than the reference approach with the constant Grüneisen parameter. The relative errors for the absorption are larger and the relative errors for the scattering are smaller when compared to the the reference approach with the correct Grüneisen parameter distribution. Further, the relative errors are smaller when compared with the reference approach with the constant Grüneisen parameter. Thus, the simulations show that simultaneous estimation of the absorption, scattering and Grüneisen parameter may be possible using the QPAT augmented with surface light measurements.

5. Conclusions

In this work, improving QPAT by combining information from surface measurements of light was investigated. The QPAT image reconstruction problem was formulated as a minimisation problem in which absorption and scattering distributions were reconstructed utilising both QPAT and surface light data. Furthermore, simultaneous estimation of the Grüneisen parameter was investigated. The approach was tested with simulations.

The results show that including surface light measurements into the QPAT image reconstruction improves both quality of the reconstructed absorption and scattering images as well as their quantitative estimates. The results are in line with other approaches in which surface light data has been utilised in QPAT [46]. Furthermore, the simulations demonstrate some of the situations in which the augmented QPAT is most useful when compared to the conventional QPAT, i.e. in distinguishing scattering, in larger domains, and when the medium consists of fine structures.

The augmented QPAT was also tested in simultaneous estimation of the optical parameters and the Grüneisen parameter. The results show that simultaneous estimation of all parameters results into better estimates for absorption and scattering than the conventional QPAT approach where the Grüneisen parameter was assumed to be constant. This indicates that it may be possible to estimate the Grüneisen parameter simultaneously with the optical parameters if QPAT is augmented with surface light measurements. However, it should be noted that the uniqueness of this minimisation problem remains to be studied.

In this work, the augmented QPAT was tested with numerical simulations in a two-dimensional setting using the diffusion approximation as the light transport model at one optical wavelength. The future work includes implementing the methodology using more realistic and numerically more challenging radiative transfer equation as well as extension of the methodology to spectral QPAT.

Funding

Academy of Finland (projects 286247 and 250215 Finnish Centre of Excellence in Inverse Problems Research), Magnus Ehrnrooth foundation, and Jane and Aatos Erkko foundation.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77, 041101 (2006). [CrossRef]

2. C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol. 54, R59–R97 (2009). [CrossRef] [PubMed]

3. L. V. Wang, ed., Photoacoustic Imaging and Spectroscopy (CRC Press, 2009). [CrossRef]

4. P. Beard, “Biomedical photoacoustic imaging,” Interface Focus 1, 602–631 (2011). [CrossRef]

5. J. Xia and L. V. Wang, “Small-animal whole-body photoacoustic tomography: a review,” Phys. Med. Biol. 61, 1380–1389 (2014).

6. L. V. Wang and J. Yao, “A practical guide to photoacoustic tomography in the life sciences,” Nature Methods 13, 627–638 (2016). [CrossRef] [PubMed]

7. J. Weber, P. C. Beard, and S. Bohndiek, “Contrast agents for molecular photoacoustic imaging,” Nature Methods 13, 639–650 (2016). [CrossRef] [PubMed]

8. J. Brunker, J. Yao, J. Laufer, and S. E. Bohndiek, “Photoacoustic imaging using genetically encoded reporters: a review,” J. Biomed. Opt. 22, 070901 (2017). [CrossRef]

9. B. Cox, J. G. Laufer, S. R. Arridge, and P. C. Beard, “Quantitative spectroscopic photoacoustic imgaging: a review,” J. Biomed. Opt. 17, 061202 (2012). [CrossRef]

10. B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am A 26, 443–455 (2009). [CrossRef]

11. D. Razansky, J. Baeten, and V. Ntziachristos, “Sensitivity of molecular target detection by multispectral optoacoustic tomography (MSOT),” Med. Phys. 36, 939–945 (2009). [CrossRef] [PubMed]

12. J. Laufer, B. Cox, E. Zhang, and P. Beard, “Quantitative determination of chromophore concentrations form 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt. 49, 1219–1233 (2010). [CrossRef] [PubMed]

13. G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in a diffusive regime,” Inv. Probl. 28, 025010 (2012). [CrossRef]

14. D. Razansky, “Multispectral optoacoustic tomography - volumetric color hearing in real time,” IEEE Sel.Top. Quantum Electron. 18, 1234–1243 (2012). [CrossRef]

15. A. V. Mamonov and K. Ren, “Quantitative photoacoustic imaging in radiative transport regime,” Comm. Math. Sci. 12, 201–234 (2014). [CrossRef]

16. A. Pulkkinen, B. T. Cox, S. R. Arridge, J. P. Kaipio, and T. Tarvainen, “A Bayesian approach to spectral quantitative photoacoustic tomography,” Inv. Probl. 30, 065012 (2014). [CrossRef]

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

18. 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, 1866–1875 (2006). [CrossRef] [PubMed]

19. 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, 2347–2356 (2008). [CrossRef]

20. T. Jetzfellner, D. Razansky, A. Rosenthal, R. Schulz, and K. H. Englmeier, “Performance of iterative optoacoustic tomography with experimental data,” Appl. Phys. Lett. 95, 013703 (2009). [CrossRef]

21. L. Yao, Y. Sun, and H. Jiang, “Quantitative photoacoustic tomography based on the radiative transfer equation,” Opt. Lett. 34, 1765–1767 (2009). [PubMed]

22. H. Gao, H. Zhao, and S. Osher, “Bregman methods in quantitative photoacoustic tomography,” UCLA CAM Report10–24 (2010).

23. R. J. Zemp, “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt. 49, 3566–3572 (2010). [CrossRef] [PubMed]

24. G. Bal and K. Ren, “Multi-source quantitative photoacoustic tomography in a diffusive regime,” Inv. Probl. 27, 075003 (2011). [CrossRef]

25. T. Tarvainen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, “Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography,” Inv. Probl. 28, 084009 (2012). [CrossRef]

26. S. Bu, Z. Liu, T. Shiina, K. Kondo, M. Yamakawa, K. Fukutani, Y. Someda, and Y. Asao, “Model-based reconstruction integrated with fluence compensation for photoacoustic tomography,” IEEE Trans. Biomed. Eng. 59, 1354–1363 (2012). [CrossRef] [PubMed]

27. X. L. Deán-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imag. 31, 1922–1928 (2012). [CrossRef]

28. T. Tarvainen, A. Pulkkinen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, “Bayesian image reconstruction in quantitative photoacoustic tomography,” IEEE Trans. Med. Imag. 32, 2287–2298 (2013). [CrossRef]

29. A. Pulkkinen, V. Kolehmainen, J. P. Kaipio, B. T. Cox, S. R. Arridge, and T. Tarvainen, “Approximate marginalization of unknown scattering in quantitative photoacoustic tomography,” Inv. Probl. Imag. 8, 811–829 (2014). [CrossRef]

30. W. Naetar and O. Scherzer, “Quantitative photoacoustic tomography with piecewise constant material parameters,” SIAM J. Imaging Sci. 7, 1755–1774 (2014). [CrossRef]

31. X. Zhang, W. Zhou, X. Zhang, and H. Gao, “Forward-backward splitting method for quantitative photoacoustic tomography,” Inv. Probl. 30, 125012 (2014). [CrossRef]

32. E. Malone, S. Powell, B. T. Cox, and S. Arridge, “Reconstruction-classification method for quantitative photoacoustic tomography,” J. Biomed. Opt. 20, 126004 (2015). [CrossRef] [PubMed]

33. A. Hannukainen, N. Hyvönen, H. Majander, and T. Tarvainen, “Efficient inclusion of total variation type priors in quantitative photoacoustic tomography,” SIAM J. Imag. Sci. 9, 1132–1153 (2016). [CrossRef]

34. M. Venugopal, P. van Es, S. Manohar, D. Roy, and R. M. Vasu, “Quantitative photoacoustic tomography by stochastic search: direct recovery of the optical absorption field,” Opt. Lett. 41, 4202–4205 (2016). [CrossRef] [PubMed]

35. R. Hochuli, S. Powell, S. Arridge, and B. Cox, “Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance,” J. Biomed. Opt. 21, 126004 (2016). [CrossRef] [PubMed]

36. P. Shao, T. Harrison, and R. J. Zemp, “Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data,” Biomed. Opt. Express 3, 3240–3248 (2012). [CrossRef] [PubMed]

37. N. Song, C. Deumié, and A. Da Silva, “Considering sources and detectors distributions for quantitative photoacoustic tomography,” Biomed. Opt. Express 5, 3960–3974 (2014). [CrossRef] [PubMed]

38. M. Haltmeier, L. Neumann, and S. Rabanser, “Single-stage reconstruction algorithm for quantitative photoacoustic tomography,” Inv. Probl. 31, 065005 (2015). [CrossRef]

39. H. Gao, J. Feng, and L. Song, “Limited-view multi-source quantitative photoacoustic tomography,” Inv. Probl. 31, 065004 (2015). [CrossRef]

40. T. Ding, K. Ren, and S. Vallélian, “A one-step reconstruction algorithm for quantitative photoacoustic imaging,” Inv. Probl. 31, 095005 (2015). [CrossRef]

41. A. Pulkkinen, B. T. Cox, S. R. Arridge, H. Goh, J. P. Kaipio, and T. Tarvainen, “Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography,” IEEE Trans. Med. Imag. 35, 2497–2508 (2016). [CrossRef]

42. B. Cox, T. Tarvainen, and S. Arridge, “Multiple illumination quantitative photoacoustic tomography using transport and diffusion models,” in “Tomography and Inverse Transport Theory (Contemporary Mathematics),”, G. Bal, D. Finch, P. Kuchment, J. Schotland, P. Stefanov, and G. Uhlmann, eds., vol. 559, pp. 1–12, (American Mathematical Society, Providence, 2011). [CrossRef]

43. X. Li, L. Xi, R. Jiang, L. Yao, and H. Jiang, “Integrated diffuse optical tomography and photoacoustic tomography: phantom validations,” Biomed. Opt. Express 2, 2348–2353 (2011). [CrossRef] [PubMed]

44. X. Li and H. Jiang, “Impact of inhomogeneous optical scattering coefficient distribution on recovery of optical absorption coefficient maps using tomographic photoacoustic data,” Phys. Med. Biol. 58, 999–1011 (2013). [CrossRef] [PubMed]

45. C. Xu, P. D. Kumavor, A. Aguirre, and Q. Zhu, “Investigation of a diffuse optical measurements-assisted quantitative photoacoustic tomographic method in reflection geometry,” J. Biomed. Opt. 17, 061213 (2012). [CrossRef] [PubMed]

46. K. Ren, H. Gao, and H. Zhao, “A hybrid reconstruction method for quantitative PAT,” SIAM J. Imag. Sci. 6, 32–55 (2013). [CrossRef]

47. T. Saratoon, T. Tarvainen, B. T. Cox, and S. R. Arridge, “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Inv. Probl. 29, 075006 (2013). [CrossRef]

48. A. Pulkkinen, B. T. Cox, S. R. Arridge, J. P. Kaipio, and T. Tarvainen, “Quantitative photoacoustic tomography using illuminations from a single direction,” J. Biomed. Opt. 20, 036015 (2015). [CrossRef] [PubMed]

49. B. A. Kaplan, J. Buchmann, S. Prohaska, and J. Laufer, “Monte-Carlo-based inversion scheme for 3D quantitative photoacoustic tomography,” in “Photons Plus Ultrasound: Imaging and Sensing 2017, Proc of SPIE,”, A. Oraevsky and L. Wang, eds., vol. 10064, pp. 100645J, (2017).

50. A. Ishimaru, Wave Propagation and Scattering in Random Media, vol. 1 (Academic Press, New York, 1978).

51. B. T. Cox and P. C. Beard, “Fast calculation of pulsed photoacoustic fields in fluids using k-space methods,” J. Acoust. Soc. Am. 117, 3616–3627 (2005). [CrossRef] [PubMed]

52. S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl. 15, R41–R93 (1999). [CrossRef]

53. S. R. Arridge and M. Schweiger, “Direct calculation of the moments of the distribution of photon time of flight in tissue with a finite-element method,” Appl. Opt. 34, 2683–2687 (1995). [CrossRef] [PubMed]

54. J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems (Springer, New York, 2005).

55. S. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20, 299–309 (1993). [CrossRef] [PubMed]

56. C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (MIT Press, 2006).

		nodes	elements
varying domain size	20 mm × 20 mm:	3975	7708
	40 mm × 40 mm:	4089	7936
	60 mm × 60 mm:	3750	7258

varying parameter distributions (20 mm × 20 mm)	2 cubes:	3673	7104
	8 cubes:	3881	7520
	18 cubes:	3717	7192
	32 cubes:	4010	7778

simultaneous estimation of the Grüneisen parameter	(10 mm × 10 mm):	3966	7690

domain size	nodes	elements
10 mm × 10 mm:	2013	3784
20 mm × 20 mm:	2489	4736
40 mm × 40 mm:	2489	4736
60 mm × 60 mm:	2218	4194

	η_{μ_a}	σ_{μ_a}	$η_{μ_{s}^{'}}$	$σ_{μ_{s}^{'}}$	η_G	σ_G	ξ
varying domain size:	0.1	0.2	2	0.5			1
varying parameter distributions:	0.1	0.2	2	0.5			1
simultaneous estimation of the Grüneisen:	0.2	0.2	1	0.25	0.3	0.1	1

		augmented QPAT		conventional QPAT

		E_{μ_a}	$E_{μ_{s}^{'}}$	E_{μ_a}	$E_{μ_{s}^{'}}$
varying domain size	20 mm × 20 mm:	4.4	18.0	8.9	25.0
	40 mm × 40 mm:	7.5	26.9	26.1	33.9
	60 mm × 60 mm:	22.4	28.6	31.2	34.2

varying parameter distributions	2 cubes:	4.2	20.3	19.1	34.7
	8 cubes:	8.7	31.5	18.5	37.5
	18 cubes:	7.7	34.9	24.5	42.1
	32 cubes:	7.4	37.5	26.4	44.2

	E_{μ_a}	$E_{μ_{s}^{'}}$	E_G
augmented QPAT:	9.49	6.09	7.04
conventional QPAT, constant Grüneisen parameter:	20.63	15.04
conventional QPAT, correct Grüneisen parameter:	4.71	26.39

Quantitative photoacoustic tomography augmented with surface light measurements

Abstract

1. Introduction

2. Methods

2.1. Forward model

2.1.1. Measuring surface light

2.2. Inverse problem of QPAT

2.2.1. QPAT augmented with surface light data

2.2.2. Simultaneous estimation of the Grüneisen parameter

3. Numerical implementations

3.1. FE-approximation of the DA

3.2. Gauss-Newton iteration

4. Results

4.1. Data simulation

4.2. Reconstructing the parameters of interest

4.3. Results

4.3.1. Varying domain size

4.3.2. Variations in the optical parameter distribution

4.3.3. Simultaneous estimation of the Grüneisen parameter

5. Conclusions

Funding

Disclosures

References and links

Cited By

Figures (4)

Tables (5)

Equations (44)

Biomedical Optics Express