Reconstruction-classification method for quantitative photoacoustic tomography

Emma Malone; Samuel Powell; Ben T. Cox; Simon Arridge

doi:10.1117/1.JBO.20.12.126004

14 December 2015 Reconstruction-classification method for quantitative photoacoustic tomography

Emma Malone, Samuel Powell, Ben T. Cox, Simon Arridge

Author Affiliations +

Journal of Biomedical Optics, Vol. 20, Issue 12, 126004 (December 2015). https://doi.org/10.1117/1.JBO.20.12.126004

Abstract

We propose a combined reconstruction-classification method for simultaneously recovering absorption and scattering in turbid media from images of absorbed optical energy. This method exploits knowledge that optical parameters are determined by a limited number of classes to iteratively improve their estimate. Numerical experiments show that the proposed approach allows for accurate recovery of absorption and scattering in two and three dimensions, and delivers superior image quality with respect to traditional reconstruction-only approaches.

1. Introduction

Photoacoustic tomography (PAT) is an emerging technique for in vivo imaging of soft biological tissue.¹ This hybrid modality uses ultrasound to detect optical contrast, combining the high resolution of acoustic methods with the spectroscopic capability of optical imaging. To generate a PA image, a short laser pulse is shined into the object, the ultrasonic waves emitted following the heating of the tissue are measured, and an image of the absorbed optical energy field is recovered. Whereas purely optical methods suffer from poor spatial resolution, acoustic waves propagate with minimal scattering, and PAT can achieve $100 μ m$ resolution at depths of several centimeters. However, PA images provide only qualitative information about the tissue, and are not directly related to tissue morphology and functionality. The principal difficulty is that the PA image is the product of both the optical absorption coefficient (which is directly related to underlying tissue composition) and the light distribution (which is not). This severely restricts the range of applications for which PAT is suitable.

Quantitative photoacoustic tomography (QPAT) aims to provide clinically valuable images of the optical absorption and scattering coefficients, or chromophore (light-absorbing molecules) concentrations from conventional PA images via an image reconstruction method.²^,³ A model of light propagation is required to relate the absorbed optical energy to the light fluence and tissue parameters. The primary challenge of QPAT is solving the nonlinear imaging problem. In particular, recovering the scattering coefficient is especially difficult due to the weak dependence of the absorbed energy density on scattering.

In this paper, we develop a method for solving the image reconstruction problem for QPAT by alternating reconstruction and segmentation steps in an automated iterative process. We introduce a probabilistic model that describes optical properties in terms of a limited number of optically distinct classes, which may correspond to tissues or chromophores. These are identified and characterized by a classification, or segmentation, algorithm. This approach allows for the use of information retrieved by the classification in the reconstruction stage and vice versa. The aim of the reconstruction is to choose solutions for which the image parameters take values close to a finite set of discrete points. The aim of the classification algorithm is to progressively improve the parametric optical model and correct for errors in the initial assumptions. Multinomial models have been employed previously in the related fields diffuse optical tomography⁴ and electrical impedance tomography.⁵ For QPAT, the main advantage is that this approach enables accurate recovery of both the absorption and scattering coefficients simultaneously.

2. Numerical Methods

2.1.

Quantitative Photoacoustic Imaging

A conventional PAT image is proportional to the absorbed optical energy

Eq. (1)

H (r) = \hat{Γ} (r) μ_{a} (r) φ [μ_{a} (r), μ_{s}^{'} (r)] r \in Ω,

where

r

is a position vector within the domain

Ω

,

μ_{a}

and

μ_{s}^{'}

are the optical absorption and reduced scattering coefficients,

φ

is the optical fluence, and

\hat{Γ}

is the Grüneisen parameter. The Grüneisen parameter represents the efficiency with which the tissue converts heat into acoustic pressure, and is often taken to be the constant

\hat{Γ} (r) = 1

\forall r \in Ω

. The fluence is dependent on the optical parameters and illumination pattern in the whole domain. The problem of recovering the optical parameters

(μ_{a}, μ_{s}^{'})

from a conventional PAT image is known as the “quantitative” problem. The optical absorption

μ_{a}

is of particular interest because it is fundamentally related to underlying tissue physiology and functionality, and encodes clinically useful information such as tissue oxygenation levels and chromophore concentrations. Conversely, the absorbed energy density

H

depends nontrivially on optical absorption and thus is not directly related to tissue morphology because it is distorted structurally and spectrally by the nonuniform light fluence.

2.2.

Diffusion Model of Light Transport

In order to recover the optical parameters $(μ_{a}, μ_{s}^{'})$ , a model of light propagation within the tissue is required. For highly scattering media and those far from boundaries and sources, a low-order spherical harmonic approximation to the “radiative transfer equation” is suitable. The “diffusion approximation” is given by⁶

Eq. (2)

[μ_{a} - \nabla \cdot κ (r) \nabla] φ (r) = q (r),

where

q (r)

is an isotropic source term and

κ = 1 / 3 μ_{s}^{'}

is the diffusion coefficient.

We set Robin boundary conditions

Eq. (3)

φ (r) + \frac{1}{2 A} κ (r) \hat{n} \cdot \nabla φ (r) = 0 r \in δ Ω,

where

A

accounts for the refractive index mismatch at the boundary.

2.3.

Minimization-Based Quantitative Photoacoustic Tomography Imaging

In this paper, we adopt a gradient-based minimization approach to image reconstruction. Typically, both $μ_{a}$ and $μ_{s}^{'}$ are unknown and need to be recovered simultaneously from the absorbed energy density. An objective function is defined, which measures the distance between the conventional PAT image $H^{m}$ and the data predicted by the model for the current estimates $H (μ_{a}, μ_{s}^{'})$

Eq. (4)

E = \frac{1}{2} \int_{Ω} {[H^{m} - H (μ_{a}, μ_{s}^{'})]}^{2} d Ω .

In order to treat the problem for a generic geometry, the finite element method is employed, whereby a weak formulation of the diffusion approximation [Eq. (2)] is considered. A discretization of the domain is defined, and the fluence and optical parameters are expressed in terms of the same piecewise-linear basis functions

u_{i} (r)

:

χ \approx \sum_{i} χ_{i} u_{i} (r)

for

χ \in {μ_{a}, μ_{s}^{'}, φ}

, where

χ_{i}

are nodal coefficients and

i = 1, \dots, N

.

We assume that the data $d^{m}$ is the absorbed energy density $H^{m}$ , projected onto a particular basis ${Ψ_{j}}$ ,

Eq. (5)

d^{m} = {d_{j}^{m}, j = 1, \dots, N}, d_{j}^{m} = \int_{Ω} H^{m} (r) Ψ_{j} (r) d Ω = ⟨ Ψ_{j}, H^{m} ⟩ .

Choices for

{Ψ_{j}}

include:

1. Point sampling $Ψ_{j} (r) = δ (r - r_{j})$ ,
2. Piecewise-linear sampling $Ψ_{j} = u_{j}$ ,
3. Sinc sampling $Ψ_{j} = sinc (| r - r_{j} |)$ .

Substituting into the objective function [Eq. (4)] leads to the discrete form of the objective function

Eq. (6)

E = \frac{1}{2} \sum_{j} {[d_{j}^{m} - ⟨ Ψ_{j}, H (μ_{a}, μ_{s}^{'}) ⟩]}^{2} = \frac{1}{2} \sum_{j} {[d_{j}^{m} - ⟨ Ψ_{j}, μ_{a} φ ⟩]}^{2} .

If a single illumination source is used and both absorption and scattering are undetermined, the problem is ill posed.² In this study, the nonuniqueness of the solution was removed by using multiple illumination patterns,⁷^–⁹ thus the objective function must be summed over the number of sources. In the following, we have omitted this sum for ease of notation. Prior information regarding the solution can be included by adding a regularization term

Eq. (7)

E = \frac{1}{2} \sum_{j} {(d_{j}^{m} - ⟨ Ψ_{j}, μ_{a} φ ⟩)}^{2} + R (μ_{a}^{'}, μ_{s}^{'}) .

In the Bayesian framework, an image is obtained by maximizing the posterior probability of the parameters, given the data

Eq. (8)

p (μ_{a}, μ_{s}^{'} | d^{m}) \propto p (d^{m} | μ_{a}, μ_{s}^{'}) p (μ_{a}, μ_{s}^{'}) .

Under this interpretation, the regularization term

R

is given by the negative log of the prior probability distribution

Eq. (9)

R (μ_{a}, μ_{s}^{'}) = - \log p (μ_{a}, μ_{s}^{'}) .

2.4.

Gradient Calculations

Cox et al.¹⁰ have shown that, for the continuous case, the gradient of Eq. (4) with respect to $μ_{a}$ at position $r^{0}$ is given by

Eq. (10)

{\frac{\partial E}{\partial μ_{a}} |}_{r^{0}} = - {φ (H^{m} - H) |}_{r^{0}} + {φ \cdot φ^{*} |}_{r^{0}},

where the “adjoint” light field

φ^{*}

is the solution to the equation

Eq. (11)

[μ_{a} - \nabla \cdot κ (r) \nabla] φ^{*} (r) = μ_{a} (H^{m} - H) .

In the following, we derive the expression for the gradient in the discrete case. The sampled forward model can be expressed as a vector $H = {H_{j}, j = 1, \dots, N}$

Eq. (12)

H_{j} = \int_{Ω} H (r) Ψ_{j} (r) d Ω = ⟨ Ψ_{j}, H ⟩ . = \sum_{i k} μ_{a i} φ_{k} \int_{Ω} Ψ_{j} (r) u_{i} (r) u_{k} (r) d Ω = φ^{T} C^{j} μ_{a},

where

C^{j}

is a sparse matrix indexed by

i

,

k

where the support of the basis functions

Ψ_{j} (r)

,

u_{i} (r)

,

u_{k} (r)

overlap. Taking the derivative of Eq. (6) with respect to

μ_{a i}

, we have

Eq. (13)

\frac{\partial E}{\partial μ_{a i}} = - \sum_{j} (\frac{\partial H_{j}}{\partial μ_{a i}}) (d_{j}^{m} - H_{j}) .

Using the expression for the absorbed energy density [Eq. (12)],

Eq. (14)

\frac{\partial H_{j}}{\partial μ_{a i}} = e_{i}^{T} C^{j} φ + μ_{a}^{T} C^{j} \frac{\partial φ}{\partial μ_{a i}},

where

e_{i}

is a vector of zeros with a single 1 in position

i

. Substituting into Eq. (13) gives

Eq. (15)

\frac{\partial E}{\partial μ_{a i}} = - \sum_{j} (e_{i}^{T} C^{j} φ + μ_{a}^{T} C^{j} \frac{\partial φ}{\partial μ_{a i}}) (d_{j}^{m} - H_{j}) .

The first term in Eq. (15) is

Eq. (16)

\sum_{j} e_{i}^{T} C^{j} φ (d_{j}^{m} - H_{j}) = \sum_{j, i, k} e_{i} C_{i k}^{j} φ_{k} (d_{j}^{m} - H_{j}) = \sum_{j, k} φ_{k} (d_{j}^{m} - H_{j}) \int_{Ω} Ψ_{j} (r) u_{i} (r) u_{k} (r) d Ω = φ^{T} E^{i} (d^{m} - H),

where

E^{i}

is given by a reordering of

C_{i k}^{j}

Eq. (17)

E_{k j}^{i} = \int_{Ω} Ψ_{j} (r) u_{i} (r) u_{k} (r) d Ω .

Note that while

C^{j}

is symmetric, in general,

E^{i}

is not.

It remains to determine $\partial φ / \partial μ_{a i}$ . The discrete form of the DA model [Eq. (2)] assumes the form¹¹

Eq. (18)

(M + K + F) φ = Q,

where

Eq. (19)

M_{j k} = \sum_{i} μ_{a i} \int_{Ω} u_{i} u_{j} u_{k} d Ω,

Eq. (20)

K_{j k} = \sum_{i} κ_{i} \int_{Ω} u_{i} \nabla u_{j} \cdot \nabla u_{k} d Ω,

Eq. (21)

F_{j k} = \sum_{i} \frac{1}{2 A} \int_{\partial Ω} u_{j} u_{k} d S,

Eq. (22)

Q_{j} = \sum_{i} q_{i} \int_{Ω} u_{i} u_{j} d Ω .

Taking the derivative of Eq. (18) with respect to the

i

’th coefficient of

μ_{a}

,

Eq. (23)

(M + K + F) \frac{\partial φ}{\partial μ_{a i}} = - V_{μ_{a}}^{i} φ,

where

Eq. (24)

V_{μ_{a}, j k}^{i} = \int_{Ω} u_{i} u_{j} u_{k} d Ω

is given by the derivative of the system matrix. We define the adjoint field

φ^{*}

as the solution to the equation

Eq. (25)

(M + K + F) φ^{*} = Q^{*},

where

Eq. (26)

Q^{*} = \sum_{j} μ_{a}^{T} C^{j} (d_{j}^{m} - H_{j})

is the adjoint source. Taking

φ^{*}

Equation (23)

- (\partial φ / \partial μ_{a i})

. Eq. (25), we obtain

Eq. (27)

\sum_{j} μ_{a}^{T} C^{j} \frac{\partial φ}{\partial μ_{a i}} (d_{j}^{m} - H_{j}) = - φ^{T} V_{μ_{a}}^{i} φ^{*} .

Substituting into Eq. (15) gives the expression for the derivative with respect to

μ_{a i}

Eq. (28)

\frac{\partial E}{\partial μ_{a i}} = φ^{T} [V_{μ_{a}}^{i} φ^{*} - E^{i} (d^{m} - H)] .

The derivative with respect to

μ_{s i}^{'}

can be derived analogously

Eq. (29)

\frac{\partial E}{\partial μ_{s i}^{'}} = - \frac{\partial κ_{i}}{\partial μ_{s i}^{'}} φ^{T} V_{μ_{s}^{'}}^{i} φ^{*},

where

Eq. (30)

V_{μ_{s}^{'}, j k}^{i} = \int_{Ω} u_{i} \nabla u_{j} \cdot \nabla u_{k} d Ω

and

(\partial κ_{i} / \partial μ_{s i}^{'}) = - 1 / 3 μ_{s i}^{' 2}

. Note that calculation of the gradient requires only two runs of the forward model. The forward problem was solved using the Toast++ software package.¹¹

Choosing point-sampling $Ψ_{j} (r) = δ (r - r_{j})$ gives simply $C^{j} = E^{i} = I$ . In this study, we chose piecewise-linear sampling $Ψ_{j} = u_{j}$ , so we had $C^{j} = E^{i} = V_{μ_{a}}^{i}$ and

Eq. (31)

\frac{\partial E}{\partial μ_{a i}} = φ^{T} V_{μ_{a}}^{i} (φ^{*} - d^{m} + H) .

3. Reconstruction-Classification Method for Quantitative Photoacoustic Tomography

A reconstruction-classification scheme is devised, which enables the recovery $μ_{a}$ and $μ_{s}^{'}$ by approaching the image reconstruction and segmentation problems simultaneously. At each reconstruction step, we minimize a regularized objective function, where the regularization term is given by a mixture model. At each classification step, the result of the previous reconstruction step is employed to update the class parameters for the multinomial model. We alternate between reconstruction and classification steps for a fixed number of iterations (Fig. 1).

3.1.

Mixture Model for $μ_{a}$ and $μ_{s}^{'}$

In this section, we introduce a probability model for $μ_{a}$ and $μ_{s}^{'}$ , which encodes prior knowledge about the optical parameters and allows us to bias the solution of the imaging problem accordingly. We assume that an array of labels $ζ_{i}$ can be determined for each node, such that

Eq. (32)

ζ_{i j} = {\begin{cases} 1 & if the i^{'} th node is assigned to the j^{'} th class, \\ 0 & otherwise . \end{cases}

The labels constitute “hidden variables” on which the image parameters are dependent. For each class

j = 1, \dots, J

, a mean vector

m_{j} = ({\bar{μ}}_{a j}, {\bar{μ}}_{s j}^{'}) \in R^{2}

is defined, and the covariance of each class is described by matrix

Σ_{j} \in R^{2 \times 2}

.

We assume that if $ζ_{i j} = 1$ , the probability distribution for $x_{i} = (μ_{a i}, μ_{s i}^{'})$ is given by a multivariate Gaussian distribution

Eq. (33)

p (x_{i} | θ_{j}) = N (m_{j}, Σ_{j}),

where

θ_{j}

indicates the set of class parameters

(m_{j}, Σ_{j})

.

The prior probability distribution of the class properties $θ_{j}$ is given by the conjugate prior to the Gaussian distribution. Prior information about the distribution of the class means or covariances can be encoded by choosing the parameters of the conjugate prior accordingly. Using a noninformative prior for the class means we have $p (m_{j}) \propto 1$ . The conjugate prior distribution for the covariance of a normal distribution is given by the normal inverse Wishart distribution (NIW):

Eq. (34)

NIW (ν_{j}, Γ_{j}) = {| Σ_{j} |}^{- (ν + d + 1) / 2} \exp [- \frac{1}{2} Tr (Γ_{j} Σ_{j}^{- 1})],

where

d

is the dimension of the domain,

ν_{j}

indicates the number of degrees of freedom, and

Γ_{j}

is a scaling matrix. If the prior is noninformative, then

ν_{j} = 0

and

Γ_{j} = 0

, and the probability distribution of the class parameters becomes

Eq. (35)

p (θ_{j}) \propto {| Σ_{j} |}^{- (d + 1) / 2},

which is known as Jeffreys prior.

The probability that the set of labels $ζ_{i} = {ζ_{i 1}, \dots, ζ_{i j}, \dots, ζ_{i J}}$ is assigned to the $i$ ’th node is given by a multinomial distribution

Eq. (36)

p (ζ_{i} | λ) = \prod_{j} λ_{j}^{ζ_{i j}},

where

λ_{j}

is the overall probability that a node is assigned to the

j

’th class. Therefore, the joint probability for

(x_{i}, ζ_{i})

is given by the product

Eq. (37)

p (x_{i}, ζ_{i} | θ, λ) = p (x_{i} | ζ_{i}, θ) p (ζ_{i} | λ) = \prod_{j} {[λ_{j} p (x_{i} | θ_{j})]}^{ζ_{i j}} .

By marginalizing over all possible values of the indicator variables

ζ_{i j}

, a “mixture of Gaussians” model for the optical parameters is obtained

Eq. (38)

p (x_{i} | θ, λ) = \int_{ζ_{i}} p (x_{i}, ζ_{i} | θ, λ) d ζ_{i} = \sum_{j} λ_{j} p (x_{i} | θ_{j}) .

Finally, for independent nodes, the prior of the image is given by

Eq. (39)

p (x | θ, λ) = \prod_{i} \sum_{j} λ_{j} p (x_{i} | θ_{j}) .

3.1.1.

Reconstruction step

The objective function takes the form of Eq. (7), where at iteration $t$ of the reconstruction-classification algorithm, the regularization is given by Eqs. (9) and (39)

Eq. (40)

R^{t} (μ_{a}, μ_{s}^{'}) = - \log p (x | θ^{t}, λ^{t}) = - \log N (\bar{x}, Σ_{\bar{x}}) = \frac{τ}{2} {∥ L_{\bar{x}} (x - \bar{x}) ∥}^{2},

where

τ

is a regularization parameter and

Eq. (41)

{\bar{x}}_{i} = {\sum_{j} ζ_{i j} \cdot m_{j} |}_{MAP (ζ)} = m_{j^{'}} \in R^{2}

is obtained by fixing the labels to the “maximum a posteriori” estimate, given the results of the previous iteration

Eq. (42)

MAP (ζ) = \arg \max_{ζ} p (ζ | x^{t - 1}, θ^{t - 1}, λ^{t - 1}),

which is calculated in the classification step (see Sec. 3.1.2). The weighting matrix

L_{\bar{x}}

is the Cholesky decomposition of

Σ_{\bar{x}}^{- 1}

, where

Σ_{\bar{x}} \in R^{2 N \times 2 N}

is a sparse matrix of which the

i

’th

2 \times 2

block along the diagonal is

Σ_{j^{'}}

if the

i

’th element belongs to the

j

’th class.

In order to sphere the solution space, that is, to render the space dimensionless, we performed a change of variables $μ_{a} \to μ_{a} / μ_{a 0}$ and $μ_{s}^{'} \to μ_{a} / μ_{s 0}^{'}$ , where $(μ_{a 0}, μ_{s 0}^{'})$ is the initial guess for the optical parameters (in this study, we initialized to the homogeneous background). Given the size of the problem, we chose a gradient-based optimization method in order to reduce memory use and computational expense.¹² The minimization was performed using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method,¹³ with a storage memory of six iterations.

3.1.2.

Classification

The purpose of the classification step is to update the multinomial model using the result of the previous reconstruction step. First, the expected values of the labels $ζ^{t + 1}$ are computed for the current class parameters $(θ^{t}, λ^{t})$ and image $x^{t} = (μ_{a}^{t}, μ_{s}^{' t})$ (E-step). Then the model parameters are updated by maximizing the posterior probability (M-step)

Eq. (43)

p (θ, λ | x^{t}) \propto p (x^{t} | θ, λ) p (θ, λ) .

E-step: The “responsibility” $r_{i j}^{t}$ is a measure of the probability that the $i$ ’th node is assigned to the $j$ ’th class. Using Bayes’ theorem and the Gaussian mixture model [Eq. (38)], we have
Eq. (44)
$p (ζ_{i j} = 1 | x_{i}^{t}, θ^{t}, λ^{t}) = \frac{p (x_{i} | ζ_{i j} = 1, θ^{t}) p (ζ_{i j} = 1)}{p (x_{i} | θ, λ)} = \frac{λ_{j}^{t} p (x_{i}^{t} | θ_{j}^{t})}{\sum_{j} λ_{j}^{t} p (x_{i}^{t} | θ_{j}^{t})} = r_{n j}^{t} .$
The expectation for the indicator values is
Eq. (45)
$E (ζ_{i j} | x_{i}^{t}, θ^{t}, λ^{t}) = \int ζ_{i j} p (ζ_{i j} = 1 | x_{i}^{t}, θ^{t}, λ^{t}) d ζ_{i j} = 0 \times p (ζ_{i j} = 0 | x_{i}^{t}, θ^{t}, λ^{t}) + 1 \times p (ζ_{i j} = 1 | x_{i}^{t}, θ^{t}, λ^{t}) = r_{i j}^{t} .$
Therefore, the MAP estimate for the labels is
Eq. (46)
$ζ_{i j}^{t + 1} = {\begin{cases} 1 & if r_{i j}^{t} is maximum \forall j, \\ 0 & otherwise, \end{cases}$
which can be used in Eq. (42).
M-step: The parameters $(θ, λ)$ are chosen in order to maximize the log posterior
Eq. (47)
$(θ^{t + 1}, λ^{t + 1}) = \arg \max_{(θ, λ)} \log p (x^{t} | θ, λ) + \log p (θ, λ) .$
Averaging over all possible values of $ζ$ gives
Eq. (48)
$\log p (x^{t} | θ, λ) + \log p (θ, λ) = \int_{ζ} \log p (x^{t}, ζ | θ, λ) d ζ + \log p (θ, λ) .$
Using “Jensen’s inequality”¹⁴ and ignoring terms that do not depend on $(θ, λ)$ , we obtain a lower bound for the log prior
Eq. (49)
$B (θ, λ) = \sum_{i} \sum_{j} r_{i j}^{t} \log [λ_{j} p (σ_{n} | θ_{j})] + \log p (λ) + \log p (θ) = \sum_{i} \sum_{j} r_{i j}^{t} [\log (λ_{j}) + \log (| Σ_{j} |) - \frac{1}{2} {(x_{i (n)} - m_{j})}^{'} Σ_{j}^{- 1} (x_{i (n)} - m_{j})] + \sum_{j} [(α_{j} - 1) \log (λ_{j}) - \frac{ν_{j} + d + 1}{2} \log | Σ_{j} |] .$
Maximizing $B (θ, λ)$ for $\sum_{j} λ_{j} = 1$ and using noninformative priors, we obtain the update rules for the model parameters
Eq. (50)
$λ_{j}^{t + 1} = \frac{\sum_{i} r_{i j}^{t}}{N},$
Eq. (51)
$m_{j}^{t + 1} = \frac{\sum_{i} r_{i j}^{t} x_{i}}{\sum_{i} r_{i j}^{t}},$
Eq. (52)
$Σ_{j}^{t + 1} = \frac{\sum_{i} r_{i j}^{t} (x_{i} - m_{j}) {(x_{i} - m_{j})}^{T} + Γ_{j}}{\sum_{i} r_{i j}^{t} + ν_{j} + d + 1} .$

3.2.

Class Means Initialization

The number of classes $J$ and the class means $m_{j}$ were initialized by automatically segmenting the result of the first reconstruction step and averaging over the segmented areas. To segment the image [e.g., see Fig. 2(a)], we looked at a binned histogram of the image of $μ_{a}$ and chose the value $μ_{a h}$ for which the number of occurrences was highest [Fig. 2(c), column 1]. We found the first node index $h$ for which the value $μ_{a h}$ occurs, and identified the corresponding scattering value $μ_{s h}^{'}$ . Having chosen a covariance matrix $Σ_{h}$ , we computed a map of the multivariate normal probability of the $(μ_{a}, μ_{s}^{'})$ images, with mean $(μ_{a h}, μ_{s h}^{'})$ [Fig. 2(c), column 2]. A suitable choice for $Σ_{h}$ is the initial covariance of the classes. Then we selected a tolerance level ${tol}_{h}$ at which to truncate the probability map, and selected all nodes with probability higher than the tolerance as belonging to the same class as node $h$ [Fig. 2(c), column 3]. We repeated this process on the remaining nodes until all nodes were classified. Thus, the number of classes was set to the number of iterations, and the average of the optical parameters over each class was used to initialize the class means [Fig. 2(b)].

Fig. 1

Reconstruction-classification algorithm outline.

Fig. 2

Class initialization example: (a) original image of $μ_{a}$ to which we apply the segmentation; (b) result of taking average image values over the segmented areas; (c) first column, histogram of occurrences of values of $μ_{a}$ in the portion of the image requiring segmentation—value with highest number of occurrences is $μ_{a h}$ (indicated by a red cross); second column, probability density function with mean $(μ_{a h}, μ_{s h}^{'})$ and covariance $Σ_{h}$ ; third column, labels identifying nodes with probability density higher than tolerance value ${tol}_{h}$ ; each row corresponds to an iteration and a distinct class, so in this case, $J = 3$ .

3.3.

Visualization of the Results

Results obtained using the reconstruction-classification method are displayed alongside scatter plots of the nodal values recovered in the two-dimensional (2-D) feature space $(μ_{a}, μ_{s}^{'})$ [e.g., see Fig. 2(c), final column in 4]. The positions of the class means $m_{j} = ({\bar{μ}}_{a j}, {\bar{μ}}_{s j}^{'})$ are identified by a cross, and the class covariances $Σ_{j}$ are represented by ellipses. These are color-coded by class, and are indicative of the clustering of image nodal values around the class means.

4. Results

4.1.

Two-Dimensional Validation and Reconstruction

We chose a numerical phantom defined on a 2-D circular mesh with 1331 nodes and radius 25 mm. Four illumination sources were placed on the boundary at angles 0, $π / 2$ , $π$ , and $3 π / 2 rad$ . In all cases, the illumination profile was a normalized Gaussian with radius (distance from the center at which the profile drops to $1 / e$ ) 6 mm. The background optical parameters were set to $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1 {mm}^{- 1}$ . Two circular perturbations of radius 6 mm were added in positions (6 mm, 10 mm) and $(- 6 mm, - 10 mm)$ [Fig. 3(a)]. The values of the perturbations were $μ_{a} = 0.02 {mm}^{- 1}$ , $μ_{s}^{'} = 1.5 {mm}^{- 1}$ and $μ_{a} = 0.03 {mm}^{- 1}$ , $μ_{s}^{'} = 1.25 {mm}^{- 1}$ , respectively. The absorbed energy field was simulated for each illumination, and 1% white Gaussian noise was added [Fig. 3(b)]. The class covariances were initialized to

Eq. (53)

Σ_{j} = (\begin{array}{l} 10^{- 6} & 0 \\ 0 & 10^{- 1} \end{array}) \forall j = 1, \dots, 3,

where the first variable was the absorption and the second was the reduced scattering. The parameters of the Jeffreys prior were set to

Γ_{j} = Σ_{j}

\forall j

,

ν (1) = 1

for the background class and

ν (2,3) = 10

for the perturbation classes. The number of classes and optical parameters was initialized using the class means initialization method (Sec. 3.2) with

{tol}_{h} = 10^{- 5}

and

Σ_{h} = Σ_{j}

[Eq. (53)], and the labels were initialized to 1 for the background class and zero for all other classes. The tolerance of the L-BFGS algorithm was set to

tol = 10^{- 11}

, and the total number of reconstruction-classification iterations was set to

MaxIt = 10

(Fig. 4). The regularization parameter

τ = 10^{- 10}

was chosen by inspection. For comparison, images were reconstructed without introducing a prior (Fig. 5); the images were reconstructed by minimizing Eq. (6) using the L-BFGS method with

tol = 10^{- 12}

.

Fig. 3

Two-dimensional (2-D) model: (a) circular mesh and (b) absorbed energy for each illumination pattern.

Fig. 4

2-D reconstruction-classification results at iteration 1 (first row), 5 (second row), and 10 (third row). Reconstructed values of $μ_{a}$ and $μ_{s}^{'}$ (first and second columns), labels recovered for perturbation classes (third and fourth columns), and scatter plot (fifth column).

Fig. 5

2-D model and reconstruction: first column, model of $μ_{s}$ and $μ_{s}^{'}$ ; second column, reconstructed values of $μ_{a}$ and $μ_{s}^{'}$ without multinomial prior; third column, reconstructed values of $μ_{a}$ and $μ_{s}^{'}$ with multinomial prior.

4.2.

Three-Dimensional Validation and Reconstruction

We chose a three-dimensional (3-D) phantom analogous to the 2-D case, defined on a cylinder with 27,084 nodes, radius 25 mm, and height 25 mm. Two spherical inclusions of radius 6 mm were placed in (6, 10, and 0 mm) and ( $- 6$ , $- 10$ , and 0 mm) [Fig. 6(a)]. Illumination sources were Gaussian in the $x y$ -plane constant in the $z$ -axis, with radius 6 mm and length 25 mm [Figs. 6(b) and 6(c)]. PAT images were simulated for four illuminations at the cardinal points, and 1% noise was added to the absorbed energy [Fig. 6(d)]. The optical, covariance, and reconstruction parameters were set to the same values used in the 2-D case. The class initialization parameters were set to ${tol}_{h} = 10^{- 7}$ and $Σ_{h} = Σ_{j}$ . Images were reconstructed by performing 10 iterations of the reconstruction-classification method (Fig. 7).

Fig. 6

Three-dimensional (3-D) model: (a) numerical phantom and perturbation locations, (b) all illumination sources, (c) cross section of optical parameters used to simulate the data for $z = 0$ , (d) cross section of absorbed energy for each illumination pattern.

Fig. 7

3-D reconstruction-classification results at iteration 1 (first row), 5 (second row), and 10 (third row). Reconstructed values of $μ_{a}$ and $μ_{s}^{'}$ (first and second columns), labels recovered for perturbation classes (third column), and scatter plot (fourth column).

5. Discussion

5.1.

Summary of Findings

We applied the proposed reconstruction-classification algorithm to a 2-D numerical phantom with three tissues, a background, and two perturbations (Fig. 3). The optical absorption was recovered reliably within a small number of iterations, and the scattering was recovered with sufficient accuracy after approximately 10 iterations (Fig. 4). We compared the optical model with images obtained by the reconstruction-classification method and by a traditional reconstruction-only (no regularization) method (Fig. 5). We found that the reconstruction-classification method delivered superior image quality, particularly with regards to the scattering parameter. We applied the reconstruction-classification algorithm to a much larger 3-D problem (Fig. 6) and observed similar results (Fig. 7) as in the 2-D case.

5.2.

Choice of Parameters

The parametric optical model and classification algorithm introduce a number of parameters that require tuning by the user. In addition to the regularization parameter, the parameters of the Jeffreys prior $Γ$ and $ν$ and the initial guess of the class variances $Σ_{j}$ must be set before performing the classification. However, their significance is fairly intuitive, and with experience of a certain type of problem, the choice of parameters becomes natural. Visualizing the class covariance matrix $Σ_{j}$ as an ellipse, changing the value of $Γ$ varies its eccentricity, and changing $ν$ varies the length of its axes. Further, given that in the first iteration the optical absorption is recovered with superior accuracy to the scattering, it is preferable to initialize the variance of the former to a smaller value than the latter, indicating greater confidence in the imaging solution.

5.3.

Initialization of the Class Means

The purpose of the means initialization scheme is to increase automation of the method so that minimum user intervention and no prior knowledge of the number of tissues or their optical properties is required. The algorithm simply performs a segmentation of the image, then takes averages over the segmented areas to initialize the class properties (Fig. 1). Alternative segmentation techniques could have been employed; however, the advantage of the proposed approach is that it directly exploits the mixture of Gaussians model to identify the tissues. Our choice to investigate a node $h$ with $μ_{a}$ belonging to the bin with a maximum number of occurrences leads to the background tissue being identified first, followed by the perturbation tissues. The choice of the node index $h$ could have been randomized so that tissues would be identified in random order. This approach is equally valid; however, we found that in cases where tissue values were close together (such as after a single reconstruction-classification iteration), it was preferable to identify the largest classes first because the mean was estimated with greater accuracy for the classes with a larger number of samples. Further, for a given image and tolerance level, our choice renders the result of the segmentation process unique and reproducible.

5.4.

Recovery of the Scattering

From the comparison with the reconstruction-only case with no regularization (Fig. 5), it is evident that the introduction of the parametric prior enables better recovery of the scattering. The inconsistency between the quality of the recovered absorption and scattering parameters in the nonregularized case is due to the weaker dependence of the latter on the absorbed energy density with respect to the former. This results in the scattering gradient being approximately an order of magnitude smaller than the absorption gradient. Although the problem can be mitigated by sphering the solution space, variations in the data due to the scattering often fall below the noise floor. In the reconstruction-classification case, typically the absorption is recovered with good accuracy within a small number of iterations. Thus, the absorption takes values very close to the class means (resulting in small clusters), and the variance along the $μ_{a}$ direction converges to a small value. Given that the regularization term is weighted by the inverse of the covariance matrix, the dependence of the absorption gradient on the data becomes weaker at each iteration, until its magnitude is comparable or smaller to that of the scattering. In the iterations that follow, the descent of the data term of the objective function is primarily due to updates to the scattering, which converges to the correct values.

5.5.

Computational Demands

Computational performance was found to be strongly dependent on the problem size. In the 2-D case with 1331 nodes (Fig. 4), the total reconstruction time (10 outer reconstruction-classification iterations) using MATLAB on a 16-processor PC with 128 GB RAM was only 77 s. In the 3-D case with 27,084 nodes (Fig. 7), the total reconstruction time increased linearly with the number of nodes and was approximately 3.7 h on the same workstation. The increase in computation time was mostly due to much longer processing times for the L-BFGS algorithm in the reconstruction step.

5.6.

Experimental Application

In experimental situations, prior information on tissue properties may be held, such as knowledge of the characteristic optical absorption and scattering spectra of chromophores of interest. These may be obtained from the literature¹⁵ or gained through tissue sample measurements. This information could be used in one of two ways. First, a library of typical chromophores could be used to initialize the class parameters instead of the proposed class means initialization method. The classification process could then perform the function of correcting for uncertainty, errors, or local variations in the real optical properties with respect to the prior information. Alternatively, it could be used to label the chromophores found by the segmentation process and identify these as certain tissues such as, e.g., “oxygenated blood” or “fat,” on the basis of the closeness of the recovered means to the characteristic properties.

5.7.

Additional Priors

In this study, we assumed independence between nodal values; however, the mixture of Gaussian models could be used in conjunction with a spatial prior. Knowledge of smoothness or sparsity properties of the solution could be employed to introduce a homogeneous spatial regularizer such as first-order Tikhonov¹⁶ or total variation.⁷^,¹⁷ Knowledge of structural information, such as that provided by an alternative imaging method or anatomical library, could be exploited by introducing a spatially varying probability map for the optical properties.

6. Conclusions

In this paper, we proposed a method for performing image reconstruction in QPAT. We introduced a parametric class model for the optical parameters and implemented a minimization-based reconstruction algorithm. We suggested an automated method by which to initialize the parameters of the class model and proposed a classification algorithm by which to progressively update and improve those parameters after each reconstruction step. We demonstrated though 2-D and 3-D numerical examples that the reconstruction-classification method allows for the simultaneous recovery of optical absorption and scattering. In particular, we found that this approach delivered superior accuracy in the recovery of the scattering with respect to traditional gradient-based reconstruction.

Acknowledgments

This work was funded by the EPSRC Doctoral Prize Fellowship EP/M506448/1.

References

1.

P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, 1 602 –631 (2011). http://dx.doi.org/10.1098/rsfs.2011.0028 Google Scholar

2.

B. Cox et al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt., 17 061202 (2012). http://dx.doi.org/10.1117/1.JBO.17.6.061202 JBOPFO 1083-3668 Google Scholar

3.

H. Gao, S. Osher and H. Zhao, “Quantitative photoacoustic tomography,” Lect. Notes Math., 2035 131 –158 (2012). http://dx.doi.org/10.1007/978-3-642-22990-9_5 Google Scholar

4.

P. Hiltunen, S. J. D. Prince and S. Arridge, “A combined reconstruction-classification method for diffuse optical tomography,” Phys. Med. Biol., 54 6457 –6476 (2009). http://dx.doi.org/10.1088/0031-9155/54/21/002 PHMBA7 0031-9155 Google Scholar

5.

E. Malone et al., “A reconstruction-classification method for multifrequency electrical impedance tomography,” IEEE Trans. Med. Imaging, 34 (7), 1486 –1497 (2015). http://dx.doi.org/10.1109/TMI.2015.2402661 ITMID4 0278-0062 Google Scholar

6.

S. Arridge, “Optical tomography in medical imaging,” Inverse Prob., 15 R41 (1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

7.

G. Bal and K. Ren, “Multiple-source quantitative photoacoustic tomography in a diffusive regime,” Inverse Prob., 27 (7), (2011). http://dx.doi.org/10.1088/0266-5611/27/7/075003 INPEEY 0266-5611 Google Scholar

8.

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). http://dx.doi.org/10.1364/AO.50.003145 APOPAI 0003-6935 Google Scholar

9.

H. Gao, J. Feng and L. Song, “Limited-view multi-source quantitative photoacoustic tomography,” Inverse Prob., 31 (6), 065004 (2015). http://dx.doi.org/10.1088/0266-5611/31/6/065004 INPEEY 0266-5611 Google Scholar

10.

B. T. Cox, S. R. Arridge and P. C. Beard, “Gradient-based quantitative photoacoustic image reconstruction for molecular imaging,” Proc. SPIE, 6437 64371T (2007). http://dx.doi.org/10.1117/12.700031 PSISDG 0277-786X Google Scholar

11.

M. Schweiger and S. Arridge, “The Toast++ software suite for forward and inverse modeling in optical tomography,” J. Biomed. Opt., 19 040801 (2014). http://dx.doi.org/10.1117/1.JBO.19.4.040801 JBOPFO 1083-3668 Google Scholar

12.

T. Saratoon et al., “3D quantitative photoacoustic tomography using the

δ

-Eddington approximation,” Proc. SPIE, 8581 85810V (2013). http://dx.doi.org/10.1117/12.2004105 PSISDG 0277-786X Google Scholar

13.

J. Nocedal and S. Wright, “Numerical optimization,” Springer Series in Operations Research and Financial Engineering, Springer-Verlag, New York (1999). Google Scholar

14.

S. Prince, Computer Vision: Models, Learning, and Inference, Cambridge University Press, Cambridge, United Kingdom (2012). Google Scholar

15.

S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol., 58 5007 –5008 (2013). http://dx.doi.org/10.1088/0031-9155/58/14/5007 PHMBA7 0031-9155 Google Scholar

16.

T. Saratoon et al., “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Inverse Prob., 29 075006 (2013). http://dx.doi.org/10.1088/0266-5611/29/7/075006 INPEEY 0266-5611 Google Scholar

17.

T. Tarvainen et al., “Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography,” Inverse Prob., 28 084009 (2012). http://dx.doi.org/10.1088/0266-5611/28/8/084009 INPEEY 0266-5611 Google Scholar

Biographies for the authors are not available.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Emma Malone, Samuel Powell, Ben T. Cox, and Simon Arridge "Reconstruction-classification method for quantitative photoacoustic tomography," Journal of Biomedical Optics 20(12), 126004 (14 December 2015). https://doi.org/10.1117/1.JBO.20.12.126004

Published: 14 December 2015

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 11 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Scattering

Tissue optics

Absorption

Photoacoustic tomography

Image segmentation

Geometrical optics

Reconstruction algorithms

1.

Introduction

2.

Numerical Methods

2.1.

Quantitative Photoacoustic Imaging

Eq. (1)

2.2.

Diffusion Model of Light Transport

Eq. (2)

Eq. (3)

2.3.

Minimization-Based Quantitative Photoacoustic Tomography Imaging

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

2.4.

Gradient Calculations

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

Eq. (22)

Eq. (23)

Eq. (24)

Eq. (25)

Eq. (26)

Eq. (27)

Eq. (28)

Eq. (29)

Eq. (30)

Eq. (31)

3.