Considering sources and detectors distributions for quantitative photoacoustic tomography

Ningning Song; Carole Deumié; Anabela Da Silva

doi:10.1364/BOE.5.003960

1. Introduction

Photoacoustic tomography (PAT) is a multiwave imaging modality that potentially combines the high optical contrast brought by the optics and high spatial resolution by ultrasound probing. The interests in using this technique for biomedical sensing are growing since 1990s, and its applications are in breast cancer detection [1,2], small animal molecular imaging [3,4], subcutaneous structures imaging [5], functional imaging [6], and other applications.

The principle consists in illuminating a sample with a local absorbing optical contrast (vasculature, presence of a tumor for example). The absorbed energy converts into heat and the temperature of the object increases, causing thermal expansion that generates an acoustic pressure wave. The energy deposited is proportional to the fluence and the local absorption coefficient of the tissues. Hence, the primary aim of PAT is to retrieve the spatial 3D distribution of the absorption coefficient of the probed tissue. This absorption coefficient, proportional to the concentration of the different chromophores composing the tissue, is an intrinsic marker of the physiology and the metabolism. In the red and near infrared, PAT has been applied to mapping hemoglobin concentrations that can reveal the presence of a cancer, through oxy- and deoxy-hemoglobin concentrations which are markers of the oxygen consumption in the tissues [7]. Glucose monitoring with PAT has also been reported [8]. However, usually the technique is applied to probing tissues at small depths (small animals, subcutaneous examination), for which the fluence can be modelled as a Beer-Lambert’s law by considering light scattering is not perturbing too much the quantification. For deep tissues imaging, because of their high level of scattering, diffusion of light has to be taken into account in the evaluation of the fluence that depends non-linearly on both the absorption and diffusion coefficients. Isolating the absorption contribution from the energy deposited is the challenge in quantitative PAT (QPAT), as the fluence does depend itself on the optical parameters (absorption and diffusion coefficients) distributions.

In a first approximation, PAT can be limited to a purely acoustic problem: from time-resolved measurements collected at the periphery, acoustic sources are localized through the reconstruction of the initial pressure distribution, from which one can obtain the absorption coefficient distribution by making the assumption of a simple photon propagation model or a uniform distribution of light. The problem is similar to seismology passive sensing: one seeks at localizing the sources that produced the collected acoustic signal. For all these approaches, a wide field illumination is preferable in order to produce a homogenized initial pressure. To go further in the quantification, different strategies have been adopted taking advantage of specific experimental protocols.

There exist basically four approaches in QPAT which consists in recovering:

i. the absorption coefficient directly. Supposing the diffusion coefficient is known, Ripoll et al. [9] recovered the absorption map with a point source accounting for instrumental factors such as the source strength, the shape of the optical pulse, and the impulse response and finite size of the transducers; Cox et al. [10] reconstructed the absorption coefficient assuming the scattering coefficient is known a priori by use of a fixed-point iterative method; Banerjee et al. [11] recovered the absorption coefficient directly from boundary pressure measurements. Very recently, Zhou et al. [12] have introduced a calibration-free method based on the direct measurement of the counts of particles through the statistical analysis of the photoacoustic signal fluctuations.
ii. the absorption coefficient and the fluence. The fluence contribution can be filtered by making use of a contrast agent of known absorption coefficient with photoacoustic signal measurements performed before and after injection [13]. Daoudi et al. [14] proposed to combine photoacoustics with acousto-optics allowing also to overcome the contribution of the fluence in the initial pressure measurement. The fluence map can be explicitly measured through Diffuse Optical Tomography [15] (DOT). By introducing prior knowledge on the spatial frequency of the two quantities, Rosenthal et al. [16] extracted both the absorption coefficient (high frequency) and the fluence (low frequency) from the photoacoustic image. As this algorithm is not based on the explicit solution the theoretical light transport equation, it does not require explicit knowledge of the illumination geometry.
iii. the absorption and diffusion coefficients simultaneously by introducing prior knowledge (spectra of known species, diffusion spectral profile) with multiple wavelength illuminations; Bal and Ren [17] have proven that, when multiple wavelength data are available, the absorption, diffusion and Grüneisen coefficients can be reconstructed simultaneously under additional minor prior assumptions.
iv. the absorption and diffusion coefficients simultaneously by introducing active probing with multiple optical illuminations. In the last decade, a lot of efforts have been produced to improve the quantification through the uncoupled reconstructions of the different physical parameters contained in the initial pressure distribution map, that is: absorption and diffusion coefficients and the Grüneisen coefficient accounting self-consistently on the acoustic properties of the tissue [18]. Tarvainen et al. [19] reconstructed both the absorption and diffusion coefficients, by considering both the radiative transport and diffusion equations as light transport models; Bal and Ren [20] recovered the absorption, diffusion, and Grüneisen coefficient mathematically when the propagation of radiation is modelled by a second-order elliptic equation; Shao et al. [21] estimated the absorption, scattering and Grüneisen distributions from reconstructed initial pressure map with multiple optical sources; In their updated work [22], they recovered the absorption and diffusion coefficients directly from measured acoustic pressure.

With this last approach, the probing becomes active and similar with what occurs in geophysics when searching for buried objects, petroleum layers, cracks... In the present work, this approach is adopted and tested under different sources and detectors geometries. The purpose is to test its robustness under potential experimental situations with single wavelength illumination, point scanning single element transducer and limited view angle examination. The paper is organized as follows: Section 2 presents the theoretical background, the results are presented and discussed in Section 3, Section 4 is the conclusion.

2. Theory

Hereafter are recalled the definitions of the physical quantities involved in the photoacoustic phenomenon and the main lines leading to the expression of the solution of the forward model. The inverse problem is also presented. It was formulated according to the method proposed by Shao et al. [22] for multiple illuminations and detections situation.

2.1 Photoacoustic phenomenon modelling in biological tissues

The amount of heat generated by tissue is proportional to the strength of the radiation. Within the conditions of thermal and stress confinements, light pulse is typically few ns, the heating function can be treated as an instantaneous function proportional to a Dirac $δ (t)$ function:

H (r, t) = μ_{a} (r) ϕ (r, t) \approx μ_{a} (r) ϕ (r) δ (t) = E (r) δ (t) .

Here

H (r, t)

is the heating function defined as the thermal energy converted at spatial position

r

and time

t

by the electromagnetic radiation per unit volume per unit time.

ϕ (r)

is the light fluence [W.cm⁻²],

μ_{a} (r)

[cm⁻¹] is the absorption coefficient. When large tissues are considered, with illumination in the wavelength range (600-1200 nm) for which light scattering mean free path length is much shorter than the absorption mean free path length, the Diffusion Approximation holds [23] and the spatial distribution of light fluence can be obtained through the resolution of the following diffusion equation:

μ_{a} (r) ϕ (r) - \nabla . (D (r) \nabla ϕ (r)) = Q (r)

D (r)

[cm] is the diffusion coefficient and

Q (r)

[W.cm⁻³] is the illumination light source density. The absorbed energy converts into heat and the temperature of the tissue increases, and thermal expansion takes place, generating an initial acoustic pressure in the medium expressed as

p_{0} (r, u) = Γ (r) μ_{a} (r) ϕ (r, u)

,

Γ

being the Grüneisen parameter, dimensionless quantity representing the efficiency of conversion of absorbed energy to pressure, and

u = {[μ_{a} D]}^{T}

. The photoacoustic pressure

p (r, u, t)

then propagates at the speed of sound in the medium

v_{s}

as:

(\frac{\partial^{2}}{\partial t^{2}} - v_{s}^{2} \nabla^{2}) p (r, u, t) = p_{0} (r, u) \frac{\partial}{\partial t} (δ (t))

which solution in free space is well known. For simplicity, to avoid the calculation of time derivatives, the measurable quantity considered in this work is the

p^{model} (r, u, t)

defined as a quantity directly proportional to the velocity potential [4]:

p^{model} (r, u, t) = v_{s}^{2} \int_{0}^{t} d t^{'} p (r, u, t^{'}) = \int d r^{'} \frac{p_{0} (r^{'}, u)}{4 π | r - r^{'} |} δ (t - \frac{| r - r^{'} |}{v_{s}})

2.2 Inverse problem formulation

From limited number M of measured pressures, at points $r_{i}$ , due to a limited number S of sources, we want to recover the maps of the optical properties $u = {[μ_{a} D]}^{T}$ . Within a perturbation approach, one can always express any optical map distribution as $u = u_{0} + δ u$ , where $u_{0}$ is the vector representing the maps of the optical properties of a reference medium, with known optical properties distribution, and $δ u$ being considered as the vector representing the perturbation to these reference maps. If $δ u < < u_{0}$ , one can express the measurements as a Taylor series expansion, truncated here at first order (Born approximation): $p_{s}^{} (r_{d}, u_{0} + δ u, t_{τ}) = p_{s}^{model} (r_{d}, u_{0}, t_{τ}) + J_{d τ}^{s} (u_{0}) δ u$ , $J_{d τ}^{s} (u_{0}) = {\partial p_{s}^{model} (r_{d}, u, t_{τ}) / \partial u |}_{u = u_{0}}$ is the Jacobian matrix. Hence, the global forward problem, for all considered sources $s \in [1, S]$ and detectors $d \in [1, M]$ , at any time $t_{τ} \in [1, T]$ , can be formulated linearly with the following matrix form:

(J^{T} J) δ u = J^{T} Δ \Leftrightarrow H δ u = b

Here

Δ = p - p^{model}

,

H

is the Hessian. The multiplication by

J^{T}

allows handling an inversion process with smaller dimensions of the matrices involved (reduced to the dimension of the chosen reconstruction mesh). The main problem remains in expressing

J

. The medium is first meshed into L voxels of volume

Δ V_{l}

. Introducing

E_{l}^{s} = μ_{a} (r_{l}) ϕ_{s} (r_{l}, u)

and

R_{d l} = | r_{d} - r_{l} |

, the model for the considered measurable quantity at time

t_{τ} \in [0, T]

, detected by point transducer

d \in [1, M]

with position

r_{d}

, due to source

s \in [1, S]

, becomes:

p_{s}^{model} (r_{d}, u, t_{τ}) = \sum_{l = 1}^{L} \frac{Γ}{4 π R_{d l}} δ (t_{τ} - \frac{R_{d l}}{v_{s}}) E_{l}^{s} Δ V_{l}

and the Jacobian matrix:

J_{d τ}^{s} (u_{0} (r_{j})) = {\partial p_{s}^{model} (r_{d}, u, t_{τ}) / \partial u |}_{u = u_{0}} = \sum_{l = 1}^{L} α_{d l, τ}^{} β_{l, s} (u_{0} (r_{j}))

α_{d l, τ}^{} = \frac{Γ Δ V_{l}}{4 π R_{d l}} δ (t_{τ} - \frac{R_{d l}}{v_{s}})

is a time evolution transfer matrix propagating the pressure created at initial time at position

r_{l}

to point detector located at

r_{d}

;

β_{l, s} (u_{0} (r_{j}))

represents the sensitivity of the energy instantaneously deposited due to the optical perturbation

δ u

located at positions

r_{j}

,

r_{j}

belonging to the reconstruction mesh:

β_{l, s} (u_{0} (r_{j})) = {[\begin{matrix} \partial E_{l}^{s} / \partial μ_{a} (r_{j}) \\ \partial E_{l}^{s} / \partial D (r_{j}) \end{matrix}]}_{u = u_{0}} = [\begin{matrix} ϕ_{s} (r_{l}, u_{0}) δ_{l j} + μ_{a 0} (r_{l}) {\frac{\partial ϕ_{s} (r_{l}, u)}{\partial μ_{a} (r_{j})} |}_{u = u_{0}} \\ μ_{a 0} (r_{l}) {\frac{\partial ϕ_{s} (r_{l}, u)}{\partial D (r_{j})} |}_{u = u_{0}} \end{matrix}]

δ_{l j}

is the Kronecker symbol. Within the perturbation approach, knowing the Green’s function

G_{0}

of the diffusion Eq. (2):

{\begin{matrix} {\frac{\partial ϕ_{s} (r_{l}, u)}{\partial μ_{a} (r_{j})} |}_{u = u_{0}} = - G_{0} (| r_{l} - r_{j} |, u_{0} (r_{j})) ϕ (| r_{j} - r_{s} |, u_{0} (r_{j})) Δ V_{j} / D_{0} \\ {\frac{\partial ϕ_{s} (r_{l}, u)}{\partial D (r_{j})} |}_{u = u_{0}} = \nabla G_{0} (| r_{l} - r_{j} |, u_{0} (r_{j})) \cdot \nabla ϕ (| r_{j} - r_{s} |, u_{0} (r_{j})) Δ V_{j} / D_{0} \end{matrix}

In order to avoid infinite values of $G_{0}$ and $ϕ$ at positions $r_{j} = r_{l}$ , the reconstruction mesh ( $r_{j}, j \in [1, J]$ , $Δ V_{j}$ pixel surface) for the values of the optical properties was chosen to be slightly different from the modelling mesh ( $r_{l}, l \in [1, L]$ ).

2.3 Computation of the Jacobian matrix

In the following case studies, the simulations were conducted in 2D. The chosen test object (Fig. 1) was basically the same as in [22]: a 6 cm × 6 cm square, in the center area of which a 2 cm × 2 cm square was selected as the reconstruction area (Fig. 1, left). The background absorption and diffusion coefficients are respectively $μ_{a 0} = 0 .1 {cm}^{-1}$ and $D_{0} = 0 .03 cm$ . The computation area contains a diffusion perturbation $δ D = 0 .003 cm$ (centered, dimensions: 0.5cm × 0.5cm) and two absorption perturbations $δ μ_{a} = 0 .01 {cm}^{-1}$ (dimensions: 0.3cm × 1.1cm and 0.3cm × 0.3cm, 0.6 cm from the center). The perturbations were chosen relatively small (10%) in order to fulfil the Born approximation. The values correspond to optical properties measured in various situations with DOT (for instance in breast cancer imaging [24]). If the perturbations are high compared to the background optical properties, non-linear reconstruction algorithms have to be used [22]. The sources are located 0.3 cm away from the reconstruction area, both the number and the length of the sources are variable. Point detectors are located 2 cm from this computation area, right at the periphery of the object.

Fig. 1 Left: Geometry of the simulated object, light sources and transducers are placed respectively 0.3 cm and 2 cm away from the object. In the reconstruction area are placed the perturbations: two absorbers ( $δ μ_{a}$ , positions: (0.6cm;1cm) and (1.6cm;1cm), dimensions: 0.3cm × 1.1cm and 0.3cm × 0.3cm) and one diffuser ( $δ D$ , position: (1cm;1cm), dimensions: 0.5cm × 0.5cm). Right: Different meshes used in the simulations: left, FEM mesh ${r_{Δ}}$ ; middle, reconstruction square mesh ${r_{j}}$ ; right: modelling square mesh ${r_{l}}$ .

Download Full Size | PDF

2.3.1 Optical energy deposition matrix: $β_{l, s} (u_{0} (r_{j}))$

The computation of this matrix is obtained through the calculation of several matrices involving the resolution of sets of diffusion Eqs. (2): $ϕ (| r_{j} - r_{s} |, u_{0} (r_{j}))$ , and $G_{0} (| r_{l} - r_{j} |, u_{0} (r_{j}))$ and their gradients in Eq. (9), and $ϕ_{s} (r_{l}, u_{0}) = ϕ (| r_{l} - r_{s} |, u_{0} (r_{j}))$ . The diffusion Eq. (2) is solved here numerically with the Finite Element Method (FEM) with Dirichlet boundary conditions at the external frontier of the domain and Neumann boundary conditions elsewhere. The optical parameters are assigned to the different subdomains, point or line sources, with unit intensity and different lengths are distributed along the above described square. The medium is meshed with triangles ${r_{Δ}}$ (Fig. 1, right) and the equations are solved for different sources geometries in order computing $ϕ (| r_{Δ} - r_{s} |, u_{0} (r_{Δ}))$ , and $G_{0} (| r_{Δ} - r_{j} |, u_{0} (r_{j}))$ and their gradients in Eq. (9). As the meshes have to be changed with the sources positions and geometry, the resulting matrices are then projected into the regular meshes ${r_{j}}$ and ${r_{l}}$ . In practice, as mentioned above, the reconstruction mesh ${r_{j}}$ was chosen such that $r_{l} \neq r_{j}$ , hence $δ_{l j} = 0$ in (8) and one never has to calculate explicitly $ϕ_{s} (r_{l}, u_{0}) = ϕ (| r_{l} - r_{s} |, u_{0} (r_{j}))$ in (8) (that is no projection into the modelling mesh ${r_{l}}$ ). The FEM mesh typically contains 4300 elements, the modelling mesh is a regular grid of 61 × 61 pixels (pixel size: 0.1 cm), and the reconstruction mesh is 21 × 21 pixels, 1 pixel shifted in horizontal and vertical with the modelling mesh. $β_{l, s} (u_{0} (r_{j}))$ is then obtained by assembling the sensitivity matrices, and its dimension is $(S \times L) \times 2 J$ .

2.3.2 Time evolution transfer matrix propagating the pressure $α_{d l, τ}^{}$

This matrix represents the acoustic propagation between each pixel of the modelling mesh and the transducer. It carries essentially the information on the distance between detectors and the different pixels of the modelling mesh $R_{d l} = | r_{d} - r_{l} |$ . It is calculated by considering the points located at the center of each pixel of the modelling mesh ${r_{l}}$ . $Δ V_{l}$ is the surface of each pixel, The Grüneisen coefficient is chosen to be a classical one in biological tissues $Γ = 0.225$ and the speed of sound $v_{s} = 1500 m {.s}^{- 1}$ . 204 time steps are considered with a time interval of 0.2 µs. The dimension of $α_{d l, τ}^{}$ is $(M \times T) \times L$ .

The dimension of the Jacobian matrix $J$ is large $(T \times M \times S) \times 2 J$ but the Hessian $H$ dimension is reduced $2 J \times 2 J$ . Vector $Δ$ representing the discrepancy between model and measurements has dimension $(T \times M \times S) \times 1$ and b is only $2 J \times 1$ . The largest matrices $J$ and $Δ$ are computed once and are not involved in the reconstruction process.

2.4 Resolution of the inverse problem

Two different reconstruction methods were tested:

- Singular value decomposition (SVD) method for which the Hessian matrix was decomposed into: $H = U \times S \times V^{T}$ . $S$ is a diagonal matrix with same dimension as $H$ containing nonnegative singular values $σ_{i}$ ordered in decreasing order, $U$ and $V$ are unitary matrices. The inversion of (5) becomes: $δ u = V S^{- 1} U^{T} b$ . A Tikhonov regularization was adopted to avoid instabilities caused by poor conditioning of $H$ : the values of the diagonal of $S^{- 1}$ are replaced by $σ_{i} / (σ_{i}^{2} + λ^{2})$ , $λ$ is the Tikhonov regularization parameter.
- The algebraic reconstruction technique (ART) [25] that searches iteratively for the solution $δ u$ such that: $δ u^{n + 1} = δ u^{n} + λ_{n} [(b - 〈 H, δ u^{n} 〉) {‖ H ‖}^{- 2}] H$ , $λ_{n}$ is a hyperparameter that can be tuned in order to regularize or to speed up the reconstruction. In all what follows, $λ_{n}$ was kept constant and equal to 1. Non negative values constraint is also introduced in the algorithm.

These reconstruction methods are tested for different noise levels. For small noise levels, SVD provides sufficiently good results but is known to be unstable and very sensitive to noise, while the ART behaves better, especially when constrained.

3. Simulation study on sources and detectors geometries

Sets of synthetic data were generated by solving the forward problem and introducing a fraction of random noise with normal distribution. Hereafter are presented the results of simulations performed on the synthetic phantom (Fig. 1) with the above-described method. The purpose of the study is to highlight the effects of different sources and detectors distributions and geometries on the reconstruction of the perturbation of the optical parameters $δ u = {[δ μ_{a} δ D]}^{T}$ , in terms of accuracy and robustness to noise. To that end, the following aspects were considered: (i) using different reconstruction methods with different noise level: this preliminary step allowed to select the proper reconstruction algorithm; (ii) number of point sources distribution: the point source illumination has been specifically studied in [22], here the results are shown as reference; (iii) point detectors distributions: the purpose is to understand which transducers distribution would be the best, for a given number of measurements; (iv) using extended sources instead of the point sources: point sources have been studied but, in practice, extended sources are easier to handle experimentally; (v) limited angle examination: this corresponds to actual experimental situations, here the transmission geometry, where sources and detectors belong to different half-spaces, is studied.

3.1 Different reconstruction methods

Four point sources and sixty detectors (15 detectors on each side of the object) are considered for this study. A first point that can be stressed is the poor conditioning of the Hessian: after performing the SVD, and examining the spectrum of the singular values $σ_{i}$ of H, the condition number obtained for this configuration is $c o n d = \max (σ_{i}) / \min (σ_{i})$ = 1.46 × 10¹⁰. Reconstructions are expected to be unstable and highly sensitive to noise.

The two different reconstruction methods, SVD (with regularization when indicated) and ART (with non-negative constraint, 1000 iterations, and $λ_{n} = 1, \forall n$ ) were tested under different noise levels (0 up to 10⁻²). Figure 2 (Left) shows the maps of the reconstructed perturbations on the absorption ( $δ μ_{a}$ ) and diffusion ( $δ D$ ) coefficients through the different reconstruction methods with increasing noise levels. On Fig. 2 (Right) are represented cross-plots extracted from these maps along a horizontal line through the middle of the reconstruction area.

Fig. 2 Comparison of different reconstruction methods at different noise levels: Left, results of the reconstructions of the perturbations on the absorption $δ μ_{a}$ and diffusion $δ D$ coefficient maps; Right: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Download Full Size | PDF

As expected, reconstructions are highly sensitive to noise, a high value of the regularization parameter has to be introduced in the SVD, showing only few elements of H out of noise may be useful for the reconstruction, even at low noise level. With the introduction of non-negative constraint, ART seems to have better potential in the recovery of the optical parameters at higher noise levels. It is interesting to notice that, though ART reconstructs the perturbations with a higher contrast. In all what follows, ART was chosen as the reconstruction algorithm.

Hence, the aim of the following studies is to improve the conditioning of H such that the reconstructions are less sensitive to noise, by optimising the data acquisition geometry.

3.2 Sources distributions

In this study, 60 point detectors were distributed around the object and the sources number was varied (2, 4, 8 and 16 sources considered). Examining the condition number of H as a function of the number of sources (Fig. 3), the reconstructions should be indeed more stable as the number of sources increases. Although the inversion is mathematically demonstrated, 2 sources are not enough to obtain uncoupled robust reconstructions. Figure 3 Left Top shows the reconstructions obtained: with 2 sources the diffusion coefficient perturbation is almost not detected, but interestingly, the absorption perturbations are well recovered. The gain obtained in the reconstructions can be quantified through the quadratic error (QE, Fig. 3 Left Bottom) between the reconstructed values $\hat{X} = {\hat{δ μ_{a}}, \hat{δ D}}$ and the target values $X = {δ μ_{a}, δ D}$ , in the reconstruction area: $Q E = \sum_{i, j} {(\hat{X_{i j}} - X_{i j})}^{2} / \sum_{i, j} X_{i j}^{2}$ . The error in the reconstruction values decreases from 1 with 2 sources to 0.4 with 4. An improvement of a factor 10 in the condition number (Fig. 3, Right Top) is obtained by adding only two more sources. However, the improvement in increasing the number of sources might be reduced after 8 sources and the gain in the reconstructions might be small even if the number of measurements (and hence the duration of the examination) is increased. This is all-the-more visible in the cross-plots shown in Fig. 3 Right Bottom: the values reconstructed become closer to the target values, the edges of the objects are sharper and the reconstruction artifacts decrease.

Fig. 3 Influence of the number of sources. Left Top: reconstructed perturbation absorption $δ μ_{a}$ and diffusion $δ D$ coefficient maps with (a) 2, (b) 4, (c) 8 and (d) 16 point sources illuminations. Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations as a function of the number of point sources. Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right). Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Download Full Size | PDF

3.3 Detectors distributions

Four sources and 60 point detectors are now considered, the number of measurements is now kept constant but the distribution of the detectors is varied. Five situations, schematized in Fig. 4, were examined: a) even and sparse distribution around the object; b) 4-sides examination but short spacing between detectors; c) same situation but with even shorter (compressed) spacing; d) 2-sides examination horizontal detection; e) 2-sides examination vertical detection.

Fig. 4 Influence of the detectors distributions. Left Top: reconstructed perturbation absorption $δ μ_{a}$ and diffusion $δ D$ coefficient maps. (a) the detectors are evenly distributed around the object (15 detectors on each side); (b) on each side, the detectors cover a length that is half the length of the object; (c) the length covered by the detectors is half of (b); (d) 30 detectors evenly distributed on each of the two vertical sides of the object; (e) 30 detectors evenly distributed on each of the two horizontal sides of the object. Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations for the five different configurations. Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right). Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Download Full Size | PDF

Examining the reconstructions (Fig. 4 Top Left), the quadratic errors (Fig. 4 Bottom Left) and the conditioning of H (Fig. 4 Top Right) obtained for these five situations shows no major improvement in the reconstructed images. Through the values of the condition number and the quadratic errors, one may conclude that distributing the detectors evenly (Fig. 4 situation (a)) may be the most favorable situation. However, reducing the field of view of the detection to the area of interest (Fig. 4 situations (b) and (c)) seems to improve the quantification while not degrading to much the localization. Figure 4 situations (d) and (e) show explicitly the influence of the positioning of the detectors in comparison to the object geometry, especially in the rectangular absorption perturbation: situation (d), when detectors are probing the object through its thinnest width, the emitted acoustic signal is better defined in time and higher in frequency (small objects producing sharper time-resolved acoustic signals), while in situation (e), the resolution is much more degraded in the direction parallel to the detectors (crossplots, Fig. 4 Right Bottom).

3.4 Point sources versus wide field illumination

QPAT has been demonstrated by using point sources illumination [22]. However, in practice, a wide field illumination is easier and safer to handle experimentally because the irradiation dose can be decreased by distributing the energy over a larger area. For this study, four line sources illuminating the four sides of the reconstruction area were considered (Fig. 5).

Fig. 5 Influence of the sources shape. Left Top: reconstructed perturbation absorption $δ μ_{a}$ and diffusion $δ D$ coefficient maps with (a) 4 points sources; 4 line sources with a length of (b) 0.2 mm; (c) 4 mm; (d) 8 mm; (e) 26 mm; (f) full field of view probing with one single source composed of lines. Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations as a function of the length of the four sources (“Integral” represents situation (f), reconstructions obtained with a single measurement but with full field of view source illumination). Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right). Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Download Full Size | PDF

Their lengths were varied from zero, corresponding to the point source situation, up to the maximum length allowing a full angle probing of the reconstruction area. Increasing the length of the source brings major improvements in the conditioning of H (Fig. 5 Right Top) and, hence, in the reconstructions (Fig. 5 Left Top). These results show that the quantification (see QE Fig. 5 Left Bottom) and localization (see crossplots Fig. 5 Right Bottom) are both definitively improved by using wider field illumination. Hence, the number of measurements required for obtaining accurate reconstructions in terms of localization and quantification can be reduced. When comparing reconstruction results shown in Fig. 3 for situation (c) and those in Fig. 5 for situations (b) to (e), with two times less measurements, the reconstructions obtained are better with 4 line sources illuminations than with 8 point sources (QEs are sensibly the same). However, the medium has still to be probed under multiple points of views: as shown for comparison on Fig. 5 situation (f), using one single wide source does not bring sufficient information to discriminate absorption and diffusion. In this last test case, the absorption coefficient is well reconstructed while the diffusion coefficient is far underestimated, showing that this illumination protocol senses essentially absorption abnormalties (see cross-plots on Fig. 5 Right Bottom).

3.5 Example of simulations for experimental situations with limited angle examination

Examining the object under all angles may be experimentally complicated or simply not feasible. Here are reported reconstructions obtained under constrained examination geometry: the medium is probed exclusively under transmission geometry, sources and detectors can rotate and probe the object on its four sides. The purpose of these calculations is to test if satisfying reconstructions can be obtained under less favorable situations corresponding nevertheless to simpler experimental protocols. Figure 6 Top Left shows the five different synthetic experiments considered here. For all of them, the four sides of the object are probed in transmission, sequentially (this corresponds to 4 measurements with an object rotating around its center of gravity with rotation angles 0°, 90°, 180° and 270°), with (a) 1 point source illumination on one side and 15 point detectors evenly distributed on the other side; from this reference situation, the illumination source was enlarged, (b) 1 line source (length: 26 mm) and 15 detectors; then the number of detectors was reduced up to a very limited number (c) same line source and 1 detector only (in this situation, sources and detectors are interchanged compared to situation (a)); (d) same as (c) but with 2 point detectors; (e) same as (c) with 3 point detectors. The condition numbers and quadratic errors obtained for the corresponding reconstructions are also reported Fig. 6.

Fig. 6 Top Left: Schemas of the experimental situations. For each situation, only 4 measurements were performed at 0°, 90°, 180° and 270°, the rotation angles corresponding to rotations around the axis located at point O (cross marker), and perpendicular to the object plan. Schemas correspond to the measurement taken at 0°: (a) 4 points sources, 15 detectors; (b) 4 line sources with length 26 mm, 15 detectors; (c) same line sources, 1 detector; (d) same line sources, 2 detectors; (e) same line sources, 3 detectors; for each situations, four synthetic measurements were considered: these configurations (0°), and object rotated by 90°, 180° and 270°. The condition numbers of the corresponding (H) are reported for each situation. Top Right: reconstructed perturbation absorption $δ μ_{a}$ and diffusion $δ D$ coefficient maps. Bottom: Quadratic errors (QE) on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations for the five different configurations.

Download Full Size | PDF

A first remark is the reconstructions obtained with point or line sources illuminations (Fig. 6 situations (a) and (b)) are improved compared to those for which the detectors are collecting signal from all the sides of the object: when the object is probed in transmission, although the condition numbers are several orders of magnitude higher, still one can get reconstructions of satisfying quality (QE≤0.4) with highly reduced number of measurements (60 measurements only). Figure 6(c) is the oversimplified situation, a single wide source and a single point detector: the absorption coefficient is satisfyingly reconstructed while the diffusion coefficient is not. However, slightly increasing the number of detectors (2 or 3 point detectors per view Fig. 6(d) and 6(e)) allows recovering reconstructions of improved quality.

This result shows the possibility of using simplified experimental PAT setups. Reducing the number of detectors points, up to 3 only here, corresponding to 12 measurements only, of course worsens the quality of the reconstructions but still allows getting better results than with complete angle detection: in comparison with Fig. 5(e), obtained with 240 measurements, the number of measurements is reduced by a factor 20. Additionally, one may think in combining the measurements in order to improve the quality of the reconstructions: the absorption coefficient seems to be fairly well reconstructed with the simplest situation (1 source, 1 detector, Fig. 6 (c)), these results could be used as prior knowledge in a second step reconstruction using 2 or 3 additional measurements (Fig. 6(d) or 6(e)).

4. Conclusion

In the present work, a multiple illumination photoacoustic tomography algorithm has been implemented according to the method suggested in [22]. The method is quantitative in the sense that it allows reconstructing both the absorption and diffusion coefficient while the conventional approaches in PAT, with a single wide field illumination give the access to the initial pressure distribution map. A description with explicit practical implementation details was provided. Different linear reconstruction algorithms were tested: the ART with non-negative constraints was shown to be less sensitive to noise than SVD or LSQR. Even if not involved in these reconstruction loops, the computation of the Jacobian matrix is still a burden, especially if one seeks to extend the approach to the non-linear case. Gradient-based methods may be used instead because more memory-efficient as they require only gradient information to approximate the Hessian matrix [17,26,27].

With this tool, a 2D simulation study on a variety of sources and detectors geometries situations was conducted. The sets of simulations were conducted through the examination of the conditioning of the Hessian matrix. The reconstructions highlighted several important results:

- Optical parameters are reconstructed with improved quality by using wide field illuminations: this study demonstrates that it is preferable to use enlarged sources instead of multiplying the number of point sources illuminations.
- The position of the detectors with respect to the perturbations to be reconstructed is crucial and deserves additional prospections. When comparing the influence of sources and detectors, the role of the illumination and detection schemes are nested and complex. The illumination scheme seems to play a major role in the quantification, while the detection in the localization with high resolved time-of-flight measurement. An important result is that probing the sample under transmittance geometry only provides better reconstructions than with full angular coverage detection.
- Different simple experimental situations were tested and good quality reconstructions were obtained with a number of measurements extremely reduced.

Experimental validations may support these findings that may influence future experimental setups designs. The present study has highlighted the importance of the data acquisition protocol. A smart combination between different types of measurements will undoubtedly improve the quality of the reconstructions. An example of such a process would be to combine one measurement taken under a full angle illumination, with a field of illumination as wide as possible, and a set of data acquired under different illumination angles. The first measurement would serve in a first reconstruction loop that would provide initialization and constraint values for a second reconstruction loop, performed with the second series of data.

Acknowledgments

We gratefully acknowledge China Scholarship Council scholarships for graduate student Ningning Song. This work was partly supported by ANR (AVENTURES-ANR-12-BLAN-BS01-0001-04).

References and links

1. S. A. Ermilov, T. Khamapirad, A. Conjusteau, M. H. Leonard, R. Lacewell, K. Mehta, T. Miller, and A. A. Oraevsky, “Laser optoacoustic imaging system for detection of breast cancer,” J. Biomed. Opt. 14(2), 024007 (2009). [CrossRef] [PubMed]

2. M. Heijblom, D. Piras, W. Xia, J. C. G. van Hespen, J. M. Klaase, F. M. van den Engh, T. G. van Leeuwen, W. Steenbergen, and S. Manohar, “Visualizing breast cancer using the Twente photoacoustic mammoscope: What do we learn from twelve new patient measurements?” Opt. Express 20(11), 11582–11597 (2012). [CrossRef] [PubMed]

3. D. Razansky, C. Vinegoni, and V. Ntziachristos, “Multispectral photoacoustic imaging of fluorochromes in small animals,” Opt. Lett. 32(19), 2891–2893 (2007). [CrossRef] [PubMed]

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

5. H. F. Zhang, K. Maslov, and L. V. Wang, “In vivo imaging of subcutaneous structures using functional photoacoustic microscopy,” Nat. Protoc. 2(4), 797–804 (2007). [CrossRef] [PubMed]

6. X. D. Wang, Y. J. Pang, G. Ku, X. Y. Xie, G. Stoica, and L. V. Wang, “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol. 21(7), 803–806 (2003). [CrossRef] [PubMed]

7. H. F. Zhang, K. Maslov, M. Sivaramakrishnan, G. Stoica, and L. V. Wang, “Imaging of hemoglobin oxygen saturation variations in single vessels in vivo using photoacoustic microscopy,” Appl. Phys. Lett. 90, 053901 (2007).

8. J. Kottmann, J. M. Rey, J. Luginbühl, E. Reichmann, and M. W. Sigrist, “Glucose sensing in human epidermis using mid-infrared photoacoustic detection,” Biomed. Opt. Express 3(4), 667–680 (2012). [CrossRef] [PubMed]

9. 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), 031912 (2005). [CrossRef] [PubMed]

10. 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]

11. 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]

12. Y. Zhou, J. Yao, K. I. Maslov, and L. V. Wang, “Calibration-free absolute quantification of particle concentration by statistical analyses of photoacoustic signals in vivo,” J. Biomed. Opt. 19(3), 037001 (2014). [CrossRef] [PubMed]

13. J. R. Rajian, P. L. Carson, and X. Wang, “Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent,” Opt. Express 17(6), 4879–4889 (2009). [CrossRef] [PubMed]

14. K. Daoudi, A. Hussain, E. Hondebrink, and W. Steenbergen, “Correcting photoacoustic signals for fluence variations using acousto-optic modulation,” Opt. Express 20(13), 14117–14129 (2012). [CrossRef] [PubMed]

15. A. Q. Bauer, R. E. Nothdurft, T. N. Erpelding, L. V. Wang, and J. P. Culver, “Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography,” J. Biomed. Opt. 16(9), 096016 (2011). [CrossRef] [PubMed]

16. A. Rosenthal, D. Razansky, and V. Ntziachristos, “Quantitative optoacoustic signal extraction using sparse signal representation,” IEEE Trans. Med. Imaging 28(12), 1997–2006 (2009). [CrossRef] [PubMed]

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

18. G. Bal and G. Uhlmann, “Inverse diffusion theory of photoacoustics,” Inverse Probl. 26, 085010 (2010).

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

20. G. Bal and K. Ren, “Multiple-source quantitative photoacoustic tomography,” Inverse Probl. 27, 075003 (2011). [CrossRef]

21. 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]

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

23. T. Vo-Dinh, Biomedical Photonics Handbook (CRC Press, 2003).

24. V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U.S.A. 97(6), 2767–2772 (2000). [CrossRef] [PubMed]

25. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, NY, 1987).

26. H. Gao, H. Zhao, and S. Osher, “Bregman methods in quantitative photoacoustic tomography,” Univ. Calif. Los Angel. UCLA Comput. Appl. Math. Rep. Vol. 10 – 42, (2010).

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

Considering sources and detectors distributions for quantitative photoacoustic tomography

Abstract

1. Introduction