Diffuse optical tomographic imaging of biological media by time-dependent parabolic SPN equations: a two-dimensional study

Jorge Bouza Domìnguez; Yves Bérubé-Lauzière

doi:10.1117/1.JBO.17.8.086012

29 August 2012 Diffuse optical tomographic imaging of biological media by time-dependent parabolic SPN equations: a two-dimensional study

Jorge Bouza Domìnguez, Yves Bérubé-Lauzière

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 8, 086012 (August 2012). https://doi.org/10.1117/1.JBO.17.8.086012

Abstract

We investigate the problem of retrieving the optical properties (absorption and scattering) of biological tissue from a set of optical measurements. A diffuse optical tomography (DOT) algorithm that incorporates constrained optimization methods is implemented. To improve image quality, the DOT algorithm exploits full time-domain data. The time-dependent parabolic simplified spherical harmonics equations (TD-pSP _N ) are used as the forward model. Time-dependent adjoint variables are resorted to in the calculation of the gradient of the objective function. Several numerical experiments for small geometric media with embedded inclusions that mimic small animal imaging are performed. In the experiments, optical coefficient values are varied in the range of realistic values for the near-infrared spectrum, including high absorption values. Single and multiparameter reconstructions are performed with the diffusion equation and higher orders of the TD-pSP _N equations. The results suggest the DOT algorithm based on the TD-pSP _N model outperforms the DE, and accurately reconstructs optical parameter distributions of biological media both spatially and quantitatively.

1. Introduction

Optical radiation in the near-infrared (NIR) spectrum is used to probe biological tissues with quite significant success in structural and functional imaging.¹^–¹² Usually, light from a distribution of sources which has propagated through a tissue is measured at the tissue’s boundary. These measurements are then used to reconstruct a map of the internal optical properties. This medical imaging modality, known as diffuse optical tomography (DOT), has undergone fast progress pertaining to models of light propagation in biological media and algorithms for solving the inverse problem.¹^,² DOT experimental systems have been simultaneously developed with measurements being made in either the continuous wave (CW),or the time-resolved (TR) regimes, with the latter comprising frequency-domain (FD) or time-domain (TD) measurements as subsets. These different measurement schemes have collected numerous data types.¹³^–¹⁵ FD and TD data are frequently used when both scattering and absorption properties of tissues are needed. On the other hand, the uniqueness of CW-based solutions to the DOT and fluorescence DOT (FDOT) inverse problems is a topic of current discussion.¹⁶^–¹⁸ It is generally believed that TD measurements can provide the richest information content about the local distribution of optical properties in tissue.¹⁹^,²⁰ Image reconstructions based on the distinct features of TD data are able to separate absorption from scatter maps, while reducing the crosstalk effect.²¹^–²⁴ Moreover, DOT and FDOT reconstruction algorithms that use full TD data show promising results in terms of improving spatial resolution and their ability to quantitatively determine optical coefficients maps.²⁵^–²⁷ Despite the former results, DOT reconstruction algorithms are limited by the accuracy of the forward model. In the presence of small geometries and high absorption values, the diffusion equation (DE) is no longer valid as a forward model.²⁸ To overcome limitations of the DE, and computational expenses of other approaches resorting to the radiative transfer equation (RTE), or the $P_{N}$ approximation thereof, low-order transport models incorporating the simplified spherical harmonics approximation ( ${SP}_{N}$ ) were derived for biomedical optics applications.²⁹^–³² Analytical solutions to ${SP}_{N}$ -based models have appeared in the literature for simple geometries with homogeneous optical properties.³³^,³⁴ Numerical solutions for complex media were also found for simulating light propagation from external sources and fluorescent probes.²⁹^–³²^,³⁵^–³⁷ These results led to the development of ${SP}_{N}$ model-based image reconstruction algorithms for solving inverse problems in intrinsic, bioluminescence, and fluorescence DOT imaging.³⁸^–⁴¹ However, as far as we know, no image reconstruction method or algorithms based on a time-dependent ${SP}_{N}$ model was ever attempted. Such an approach could benefit from both the accuracy of the ${SP}_{N}$ models and the richness of TD data. The subject of this paper is to provide such an approach, and assess the accuracy of the TD ${SP}_{N}$ equations as a forward model (see below) for recovering the internal optical properties of biological media.

With that goal, an image reconstruction algorithm for solving the inverse DOT imaging problem is implemented. As the forward model, which is reviewed in Sec. 2, a discretized formulation of the time-dependent parabolic ${SP}_{N}$ equations (the TD- $p {SP}_{N}$ model) is used, which was derived and studied in depth in a previous work.³². The spatial discretization of the model is performed through the finite element method (FEM), and the temporal discretization is carried out with an implicit finite difference scheme. In this way, it is possible to deal with complex geometries and heterogeneous distributions of the optical coefficients. Following the latest trends in model-based image reconstruction algorithms, which resort to iterative optimization techniques,⁴²^–⁴⁶ the inverse problem is cast here as a partial differential equation (PDE)-constrained optimization scheme. This will be described in Sec. 3. To aid in reducing the ill-posedness and improve the stability of the algorithm, the inverse problem solution is further constrained by imposing bounds on the optical coefficient values. The resulting optimization problem possesses a high number of unknown variables, and has to deal with a large amount of data owing to the use of TD data. To reduce complexity, a nested analysis and design (NAND) method is employed.⁴⁷ Further, the problem is framed in a sequential quadratic programming (SQP) algorithm.⁴⁸ By the introduction of adjoint variables, an iterative numerical scheme for calculating the gradient values of the objective function (OF) at each time step is provided, with the OF measuring the discrepancy between the predicted and actual measurements. The former steps (i.e., the use of NAND, SQP, and adjoint variables) taken in solving the inverse problem reduce the computation time considerably.

To examine the DOT algorithm’s performances and the accuracy of the TD- $p {SP}_{N}$ model in the reconstructions, a series of numerical experiments resembling small animal intrinsic imaging are carried out, as described in Sec. 4. Single parameter reconstructions are performed for the case of a homogeneous medium with an embedded absorptive inclusion, for which increasing values of the absorption coefficient are employed. High absorption values which simulate practical situations encountered in DOT and FDOT are considered. The image reconstructions for ${SP}_{N}$ orders $N > 1$ yield better quantitative results than the DE-based reconstructions. Particularly, the reconstructed images using the TD- $p {SP}_{3}$ approximation provide the best results. Next, a scattering inclusion is added in the previous situations, and the simultaneous reconstruction of the absorption and the diffusion coefficient maps is performed. As before, image reconstructions for $N > 1$ better estimate the optical properties of the inclusions than the DE (TD- $p {SP}_{1}$ ).

As will be seen, the results obtained support the TD- $p {SP}_{N}$ -based DOT algorithm presented herein as an accurate approach for obtaining spatial distributions of internal optical properties of biological media.

2. Forward Problem

In Ref. 32, a new time-dependent model based on the ${SP}_{N}$ approximation for describing light propagation in scattering and absorbing media was described and validated. The model is represented by a set of $(N + 1) / 2$ -coupled parabolic partial differential equations (PDEs) referred to as the time-domain parabolic ${SP}_{N} (TD − p {SP}_{N})$ equations. The TD- $p {SP}_{N}$ equations can be written in matrix notation in terms of the vector of composite moment functions $Φ$ as³²

Eq. (1)

[C + \frac{η}{c} \frac{\partial}{\partial t} T] Φ (r, t) + D_{r} Φ (r, t) = Q (r, t),

where

η

is the refractive index of the medium, and

c

is the speed of light in a vacuum. The explicit expression for the matrices

C

,

T

, the matrix operator

D_{r}

, and the source vector

Q

are reproduced in Appendix A. FEM spatial discretization and an implicit finite differencing method (FDM) for the time variable lead to the following FEM-FDM numerical representation of the TD-

p {SP}_{N}

equations³²

Eq. (2)

\tilde{W} {\tilde{Φ}}^{(m)} = (\frac{η}{c}) \tilde{T} {\tilde{Φ}}^{(m - 1)} + {\tilde{ϒ}}^{(m)},

where

Eq. (3)

\tilde{W} = Δ t (\tilde{K} + \tilde{M} + \tilde{Π}) + (\frac{η}{c}) \tilde{T},

and

Eq. (4)

{\tilde{ϒ}}^{(m)} = Δ t [{\tilde{F}}^{(m)} + {\tilde{Γ}}^{(m)}] .

In Eqs. (2)–(4),

{\tilde{Φ}}^{(m)}

represents the vector of composite moments of the radiance at time

t_{m} = m Δ t

, where

Δ t

is the time-step employed in the finite difference scheme. The remaining terms appearing in Eqs. (3) and (4) are also described in Appendix A, along with their dependencies on the absorption

μ_{a}

, scattering

μ_{s}

, and anisotropy

g

coefficients. As the terms given in Eqs. (3) and (4) are critical in describing the solution to the inverse problem, Appendix A and Ref. 32 should be consulted as necessary. Next, the notation associated with Eq. (2) is discussed.

2.1.

Notation

The ${SP}_{N}$ equations come naturally with odd orders $N$ .²⁹^,³² To simplify the notation in the sequel, the constant $K = (N + 1) / 2$ is defined. The following notation will also be used. First, the quantity $Φ (r, t)$ considered in the TD- $p {SP}_{N}$ model is given by a vector of space- and time-dependent functions as follows

Eq. (5)

Φ (r, t) = {[\begin{matrix} φ_{1} (r, t) & φ_{2} (r, t) & \dots & φ_{K} (r, t) \end{matrix}]}^{T},

where

{[]}^{T}

means the transpose; the functions

φ_{k} (r, t)

k = 1, \dots, K

are called the composite moments, and they constitute what are referred to as the state variables when image reconstruction is considered later. When the time variable is discretized at timesteps

t_{m} = m Δ t

,

m = 1, \dots, M - 1

, the following quantities are defined

Eq. (6)

Φ^{(m)} (r) \equiv Φ (r, t_{m}) .

Using the FEM, the solutions

Φ (r, t)

to the TD-

p {SP}_{N}

equations are approximated over a mesh of elements (e.g., triangles) by a piecewise polynomial and continuous function

Φ^{h} (r, t)

. Let

L

be the number of nodes in the mesh (e.g., the different vertices of the mesh triangles). Then, the solution may be written as

Eq. (7)

Φ^{h} (r, t) = \sum_{i = 1}^{L} Φ_{i} (t) u_{i} (r), u_{i} (r) \in Ω^{h},

where

Ω^{h}

is a finite-dimensional subspace spanned by the polynomial basis functions

u_{i} (r)

,

i = 1, \dots, L

, and where

h

denotes the characteristic size of the elements. Hence, the problem is reduced to finding the nodal values

\tilde{Φ} (t) = {Φ_{i} (t)}_{i = 1, \dots, L}

from which the solution can be obtained everywhere through the interpolation rule given in Eq. (7). The simplest choice, and the one made here, is a piecewise linear basis over a mesh of triangles in which the nodal values of the basis functions are defined by

u_{i} (r_{j}) = δ_{i, j}

,

i, j = 1, \dots, L

. Let

i = 1, \dots, L

be an index over the set of all nodes in the FEM mesh of triangles that will be denoted by

D = {r_{i}}_{i = 1, \dots L}

, with

r

denoting the position of a point in space. The subset of interior nodes in

D

will be denoted by

D_{I}

, and

L_{I}

will denote the number of interior nodes. The remaining nodes form the set of boundary nodes denoted by

D_{B}

containing a total of

L_{B}

nodes, with, of course,

L_{I} + L_{B} = L

.

The discretized version of $Φ (r, t)$ is

{\tilde{Φ}}^{(m)} = {[{\tilde{φ}}_{1} (r_{1}, t_{m}) {\tilde{φ}}_{1} (r_{2}, t_{m}) \dots {\tilde{φ}}_{1} (r_{L}, t_{m}) {\tilde{φ}}_{2} (r_{1}, t_{m}) \dots {\tilde{φ}}_{2} (r_{L}, t_{m}) \dots {\tilde{φ}}_{K} (r_{1}, t_{m}) \dots {\tilde{φ}}_{K} (r_{L}, t_{m})]}^{T},

or in more abbreviated form

Eq. (8)

{\tilde{Φ}}^{(m)} = {[\begin{matrix} {\tilde{φ}}_{1}^{(m)} & {\tilde{φ}}_{2}^{(m)} & \dots & {\tilde{φ}}_{K}^{(m)} \end{matrix}]}^{T} .

Light source configurations are indexed with the subscript

s

. A source configuration is used for illuminating the medium and making a set of measurements. Such a set of measurements for a given source configuration is what will be called a “tomographic projection”. A source configuration may, for instance, consist of a set of different sources simultaneously illuminating the object. For a given source configuration

s

, the associated solution of the discretized TD-

p {SP}_{N}

equations at a given time

m

will be written as

{\tilde{Φ}}_{s}^{(m)}

. Hereon,

{{\tilde{Φ}}_{s}^{(m)}}

will denote the set containing the solution at all time steps for all source configurations.

Note that the discretized TD- $p {SP}_{N}$ equations are solved for each $s = 1, \dots, S$ (i.e., for each source configuration illuminating the object) in order to predict the measurements made for that configuration. This means that the forward problem needs to be solved as many times as there are source configurations.

In the following, $μ$ will denote the vector containing the set of optical coefficients evaluated at the nodes of the FEM mesh and to be retrieved through the solution of the inverse problem. Each optical coefficient is expanded along the same type of basis functions of the FEM expansion, see Eq. (7). In accordance with the previous notation, if for example the optical coefficients $μ_{a}$ and $D$ are considered, as will be done later, where $D$ is the diffusion coefficient given by $D = 1 / [3 (μ_{a} + (1 - g) μ_{s})]$ , then

Eq. (9)

μ = {[μ_{a} (r_{1}) μ_{a} (r_{2}) \dots μ_{a} (r_{L}) D (r_{1}) D (r_{2}) \dots D (r_{L})]}^{T}

represents the nodal values of the optical coefficients.

2.2.

Measurements Modeling

The type of measurements depends on the experimental conditions, and can take diverse forms.²² Here, measurements are taken as the full time-dependence of the outgoing normal component of the radiant current density vector $J_{n, s, d}^{(m)}$ (or excitance) at a finite collection of detector positions $r_{d}$ , $d = 1, \dots, D$ for each source configuration $s$ at $m$ different times. The values of $J_{n, s, d}^{(m)}$ can be calculated using ${\tilde{Φ}}_{s}^{(m)}$ and the optical coefficients evaluated at the detector positions via

Eq. (10)

J_{n, s, d}^{(m)} = [j_{1} - j_{2} {(B)}^{- 1} A] {\tilde{Φ}}_{s}^{(m)} (r_{d}, t) = V_{d} {\tilde{Φ}}_{s}^{(m)} (r_{d}, t), s = 1, \dots, S, d = 1, \dots, D, m = 1, \dots, M,

where the expressions for vectors

j_{1}

and

j_{2}

are given in Appendix B, and

A

and

B

are matrices appearing in the boundary conditions,²⁹^,³² see Appendix A. The elements of both matrices depend on the values of the optical coefficients of the medium at the boundary nodes. The vector

V_{d}

is a compact representation of the “measurement operator”

j_{1} - j_{2} {(B)}^{- 1} A

. Note that if the outward flux is collected in a detection area, and that if geometrical factors and the numerical aperture of a real detector are taken into account (see Sec. 7-6 in Ref. 49), then other types of measurements can be expressed as a weighted sum of the values of

J_{n, s, d}^{(m)}

.

3. Inverse Problem

DOT imaging is an inverse problem that can be considered as a parameter-estimation task. It involves finding parameter values, which are optimal in that they minimize an objective—or functional—that measures the misfit between a given set of measurements and predictions thereof. A judicious choice of functional improves the quality of the final image and accelerates the convergence of the reconstruction algorithm, even in the presence of measurement noise. In addition, a suitable functional reduces the effect of cross variations on the optical parameter subsets, an effect known as crosstalk.²² Several functionals and selection of measurement types are proposed in the literature to eliminate functional valley structure, and achieve less eccentricity in the functional characteristic surface.²²^,⁵⁰ Even so, simulated and experimental reconstructions still exhibit poor spatial resolution and quantitativeness.¹⁹^,²¹^,²² Some authors recommend exploiting the full profiles of TD measurements in nonlinear reconstruction schemes.²⁵^,²⁶ The main idea behind the mentioned approaches is that the amplitude of the TD profile is influenced by both absorption and scattering coefficients (the influence of the absorption coefficient is stronger in these respects), while the shape of the TD profileis more dependent on the scattering coefficient values.¹⁹ Still, exploiting TDdata in image reconstruction algorithms and optimal functional selection is a subject of active research.⁵⁰^–⁵²

3.1.

Optimization Problem

Following the previous discussion, a least-squares estimation is cast to the inverse problem using an objective function resembling that proposed in Ref. 25. Then, the minimum of the difference between the TD measurement data $M_{s, d}^{(m)}$ , $s = 1, \dots, S$ , $d = 1, \dots, D$ , $m = 1, \dots, M$ and the predictions $J_{n, s, d}^{(m)}$ obtained from the TD- $p {SP}_{N}$ forward model is sought. Thus, the following objective function (OF) is considered

Eq. (11)

f = \frac{1}{2} \sum_{s = 1}^{S} \sum_{d = 1}^{D} \sum_{m = 1}^{M} {[\frac{M_{s, d}^{(m)} - J_{n, s, d}^{(m)}}{σ_{s, d}^{(m)}}]}^{2},

where the magnitudes

σ_{s, d}^{(m)}

quantify the uncertainty in measurements, which are mainly influenced by photon noise. Because of the Poisson distribution nature of the photon noise present in TD data (shot noise is considered as the main contribution here); the uncertainty is assumed to be proportional to the square root of the measured signal,⁵¹ i.e.,

σ_{s, d}^{(m)} \sim \sqrt{M_{s, d}^{(m)}}

. To gain clarity for the sequel, the OF is re-written directly in terms of

{\bar{Φ}}_{s}^{(m)}

Eq. (12)

f = \frac{1}{2} \sum_{s = 1}^{S} \sum_{d = 1}^{D} \sum_{m = 1}^{M} \frac{{[M_{s, d}^{(m)} - \sum_{k = 1}^{K} V_{k, d} P_{k, d} {\bar{Φ}}_{s}^{(m)}]}^{2}}{M_{s, d}^{(m)}},

where the measurement operator

V_{d}

appearing in Eq. (10) was decomposed into

K

components (

k = 1, \dots K

). The projection operator (or vector)

P_{k, i}

is represented by a row vector of length

K

times

L

which has zero components except for a value equal to 1 at position

(k - 1) L + i

. Thus, for any

m

, [refer back to Eq. (8)]

Eq. (13)

P_{k, d} {\tilde{Φ}}_{s}^{(m)} = {\tilde{φ}}_{k, s, d}^{(m)} .

Now, the inverse problem can be posed as finding

Eq. (14)

μ = \underset{{μ^{h}}}{\arg \min} \frac{1}{2} \sum_{s = 1}^{S} \sum_{d = 1}^{D} \sum_{m = 1}^{M} \frac{{(M_{s, d}^{(m)} - \sum_{k = 1}^{K} V_{k, d} P_{k, d} {\tilde{Φ}}_{s}^{(m)})}^{2}}{M_{s, d}^{(m)}}, subject to {\begin{matrix} \tilde{W} {\tilde{Φ}}_{s}^{(m)} = (\frac{η}{c}) \tilde{T} {\tilde{Φ}}_{s}^{(m - 1)} + {\tilde{ϒ}}_{s}^{(m)}, \\ μ_{low} \leq μ \leq μ_{up}, \end{matrix} \begin{matrix} m = 1, \dots, M, s = 1, \dots, S, \end{matrix}

where

\arg \min

stands for argument of the minimum. The vectors

μ_{low}

and

μ_{up}

are lower and upper bounds over the set of optical coefficients

μ

.

The inverse problem posed in Eq. (14) is a constrained optimization problem, specifically a PDE-constrained optimization problem. In general, PDE-constrained optimization refers to the optimization of systems whose states are governed by PDEs. The forward problem consists insolvingthe PDEs for the state variables (e.g., temperature, electric field, or in the present case, the optical field ${{\tilde{Φ}}_{s}^{(m)}}$ ), given appropriate data (e.g., geometry, boundary conditions, initial conditions, source functions, and in DOT, the set of optical coefficients $μ$ ). The optimization problem seeks to determine some of these data, called the decision or design variables, given performance goals in the form of an OF [first line of Eq. (14)], and possibly inequality or equality constraints on the behaviour and parameters of the system. Since the system’s behavior is modeled by PDEs, the latter appear as equality constraints in the optimization problem. The PDE constraints are usually referred as the state equations. There is an extra feature in the present problem that increases its complexity. If the state equations are evolutionary in nature (that is, dependent on time), the optimality conditions (necessary conditions for the existence of a minimum) for Eq. (14) constitutes a boundary-value problem in the space–time domain. For this reason, the optimization problem considered here is significantly more difficult to solve than steady-state problems. Posing the inverse problem as a PDE-constrained optimization problem is a relatively new approach in biomedical optics⁴²^–⁴⁶ that benefits from considerable experience gained in the area of large-scale optimization.⁴⁷^,⁴⁸^,⁵³^,⁵⁴

3.2.

Nested Analysis and Design

The total number of variables (state and design variables) in the present case dealing with TD data considerably surpasses those in CW and FD inverse problems. The optimization problem becomes more complex because of the time-dependence of the forward model. The reconstruction algorithm has to deal with a large amount of variables and involves sparse matrices whose size and manipulation are memory demanding, especially for a fine mesh and small time-steps. Typically, a TD measurement can contain up to a few thousands of points. For the present purposes, on the order of a hundred points on each TD curve are used, which makes the number of state variables around $100 \times K times$ that of design variables. Motivated by these considerations, it was decided to employ a “nested analysis and design” (NAND) method to solve the inverse problem. The NAND method, also called the “black-box” approach, reduces the number of unknown parameters in a constrained optimization problem⁴⁷ by considering the state variables (in the present case, ${{\tilde{Φ}}_{s}^{(m)}}$ ) as implicit functions of the design variables $μ$ . Thus, the states ${{\tilde{Φ}}_{s}^{(m)}}$ are obtained by explicitly solvingthe PDE [given by Eq. (2)], and substituted in the OF, leading to a reduction of the dimensionality of the problem. In mathematical terms

Eq. (15)

{\tilde{Φ}}_{s}^{(m)} = {\tilde{Φ}}_{s}^{(m)} (μ) \Rightarrow f ({\tilde{Φ}}_{s}^{(m)}, μ) = f ({\tilde{Φ}}_{s}^{(m)} (μ), μ) = f (μ) .

Consequently, the problem is reduced to an optimization problem which is unconstrained with respect to the PDE, but constrained with respect to the optical properties to be recovered because of the lower and upper bounds imposed thereon. These upper and lower bound constraints are treated in the algorithm as linear constraints with an associated

l_{I}

merit function; see below. Descriptions of the algorithm and implementation details are provided next.

3.3.

Sequential Quadratic Programming

To find the solution, the optimization problem posed in Eq. (14) is cast into a sequential quadratic programming (SQP) algorithm. SQP basically consists in locally approximating the objective function at each iteration step by a quadratic function. This significantly simplifies the optimization problem, because a simpler quadratic function needs to be considered, rather than a nonlinear one.⁴⁸ SQP is widely used in large-scale optimization for solving nonlinear constrained optimization problems because of the reduction in the computational cost and its fast convergence. Details about the implementation of this kind of schemes can be found elsewhere. For example, see Refs. 48, 53, and 54 for extensive documentation about the topic, and Refs. 44 45.–46 and 54 for applications of inverse problems in electromagnetism, DOT, and fluorescence DOT. A description of the SQP method used in the context of the present DOT algorithm will now be made.

For the line-search at every iterate $μ_{k}$ of the optical properties at step $k$ during the optimization process, the search direction $d_{k}$ in the SQP scheme is obtained by solving the quadratic subproblem⁴⁸

Eq. (16)

\min_{d \in R^{n}} \frac{1}{2} d_{k}^{T} H_{k} d_{k} + g_{k} d_{k}, subject to μ_{low} \leq μ \leq μ_{up},

where

g_{k}

denotes the gradient of the OF at step

k

, which is obtained through the forward model (see next subsection), and

H_{k}

denotes the Hessian (with respect to

μ_{k}

) of the Lagrangian function

L (μ, λ) = f (μ) - β^{T} c (μ)

[

β

is the vector of Lagrange multipliers determined by the Karush-Khun-Tucker (KKT) conditions]. Here, the bound constraints

μ_{low} \leq μ \leq μ_{up}

are rewritten as linear constraints represented by

c (μ) \geq 0

. The updated vector is given by

Eq. (17)

μ_{k + 1} = μ_{k} + α_{k} d_{k},

where

α_{k}

is a step-length parameter that provides sufficient reduction in the

l_{1}

merit function

Eq. (18)

ϕ_{ρ} = f (μ) + {ρ ‖ v (μ) ‖}_{1} .

Here,

ρ

is the penalty parameter defined to enforce a descent of

φ_{ρ}

, and

v (μ)

measures the violation of the constraints; the function

v

is defined as

Eq. (19)

{‖ ν (μ) ‖}_{1} = \sum_{n = 1}^{N_{v}} {| \max (0, μ_{n} - μ_{n, up}) | + | \max (0, μ_{n, low} - μ_{n}) |},

where

N_{v}

is the number of elements in the vector

μ

. Usually in such SQP schemes, the Hessian is updated by the damped Broyden–Fletcher–Goldfarb–Shanno (BFGS) method to avoid computing second derivatives,⁴⁸ and this method was used here. For a representation of the general algorithm and the SQP scheme, see the flowchart appearing in Fig. 1.

Fig. 1

Flowchart of the general algorithm and sequential quadratic programming (SQP) scheme.

In the following, some focus is put on the necessary calculation of the gradient of the OF. The technique based on adjoint variables avoids using finite difference schemes of the OF and associated inherent errors and extremely high computational cost.

3.4.

Gradient of the Objective Function

Calculating the gradient of the OF is undoubtedly the most difficult part of an inverse problem. This will now be tackled for the specific problem considered herein. Here, the gradient of the OF $g$ can be expanded as a sum of independent terms for each source configuration: $g = \sum_{s = 1}^{S} g_{s}$ . Considering the implicit dependence of ${\tilde{Φ}}_{s}^{(m)}$ on $μ$ , the $j$ ’th component of the gradient $g_{s, j}$ can be written as

Eq. (20)

g_{s, j} = \sum_{m = 1}^{M} \frac{d f}{d {\tilde{Φ}}_{s}^{(m)}} \frac{\partial {\tilde{Φ}}_{s}^{(m)}}{\partial μ_{j}},

where

d f / d {\tilde{Φ}}_{s}^{(m)}

denotes the derivative of the OF with respect to

{\tilde{Φ}}_{s}^{(m)}

. The explicit variation of

{\tilde{Φ}}_{s}^{(m)}

with

μ

can be obtained from Eq. (2) as

Eq. (21)

\frac{\partial {\tilde{Φ}}_{s}^{(m)}}{\partial μ_{j}} = - {\tilde{W}}^{- 1} (\frac{\partial \tilde{W}}{\partial μ_{j}}) {\tilde{Φ}}_{s}^{(m)} .

Note that in obtaining Eq. (21),

{\tilde{Φ}}_{s}^{(m - 1)}

, which appears in Eq. (2), was not differentiated with respect to

μ_{j}

, because, as shown below, backwards differentiation will be resorted to, and in this framework,

{\tilde{Φ}}_{s}^{(m - 1)}

is considered dependent on

{\tilde{Φ}}_{s}^{(m)}

. Next, an expression for

d f / d {\tilde{Φ}}_{s}^{(m)}

is derived. Note that because of Eq. (2), a simple explicit derivative of the OF with respect to

{\tilde{Φ}}_{s}^{(m)}

does not consider the dependence of

{\tilde{Φ}}_{s}^{(m)}

on previous time values of the field which comes from the FDM scheme. To resolve that issue, the core idea of adjoint differentiation (or reverse differentiation) expounded in Ref. 55 for a time-dependent DE-based inverse problem is followed. The idea is toapply the chain rule of differentiation backwards from time

m + 1

to

m

to the OF to yield

Eq. (22)

\frac{d f}{d {\tilde{Φ}}_{s}^{(m)}} = \frac{d f}{d {\tilde{Φ}}_{s}^{(m + 1)}} \frac{d {\tilde{Φ}}_{s}^{(m + 1)}}{d {\tilde{Φ}}_{s}^{(m)}} + \frac{\partial f}{\partial {\tilde{Φ}}_{s}^{(m)}},

with

Eq. (23)

\frac{d f}{d {\tilde{Φ}}_{s}^{(M)}} = \frac{\partial f}{\partial {\tilde{Φ}}_{s}^{(M)}} .

Equations (22) and (23) now take into account the time-dependence of

{\tilde{Φ}}_{s}^{(m)}

. From Eqs. (2) and, the following terms can be calculated

Eq. (24)

\frac{d {\tilde{Φ}}_{s}^{(m + 1)}}{d {\tilde{Φ}}_{s}^{(m)}} = {\tilde{W}}^{- 1} (\frac{η}{c}) \tilde{T},

Eq. (25)

\frac{\partial f}{\partial {\tilde{Φ}}_{s}^{(m)}} = \sum_{d = 1}^{D} \frac{(\sum_{k = 1}^{K} V_{k, d} P_{k, d} {\tilde{Φ}}_{s}^{(m)} - M_{s, d}^{(m)})}{M_{s, d}^{(m)}} \sum_{k = 1}^{K} V_{k, d} P_{k, d} .

Substituting Eq. (21) into Eq. (20), and Eq. (24) into Eq. (22), leads to

Eq. (26)

g_{s, j} = - \sum_{m = 1}^{M} \frac{d f}{d {\tilde{Φ}}_{s}^{(m)}} {\tilde{W}}^{- 1} (\frac{\partial \tilde{W}}{\partial μ_{j}}) {\tilde{Φ}}_{s}^{(m)}, \frac{d f}{d {\tilde{Φ}}_{s}^{(m)}} = \frac{d f}{d {\tilde{Φ}}_{s}^{(m + 1)}} {\tilde{W}}^{- 1} (\frac{η}{c}) \tilde{T} + \frac{\partial f}{\partial {\tilde{Φ}}_{s}^{(m)}} .

Equation (26) suggests the introduction of the adjoint variables

λ_{s}^{(m)}

as auxiliary variables

Eq. (27)

{\tilde{W}}^{T} λ_{s}^{(m)} = - {(\frac{d f}{d {\tilde{Φ}}_{s}^{(m)}})}^{T},

which allows rewriting Eq. (26) as

Eq. (28)

g_{s, j} = \sum_{m = 1}^{M} {(λ_{s}^{(m)})}^{T} (\frac{\partial \tilde{W}}{\partial μ_{j}}) {\tilde{Φ}}_{s}^{(m)},

Eq. (29)

\frac{d f}{d {\tilde{Φ}}_{s}^{(m)}} = {- (λ_{s}^{(m + 1)})}^{T} (\frac{η}{c}) \tilde{T} + \frac{\partial f}{\partial {\tilde{Φ}}_{s}^{(m)}} .

Introducing the adjoint variables

λ_{s}^{(m)}

reduces the number of computations to be carried out because it avoids the need to compute the partial derivatives of the field

{\tilde{Φ}}_{s}^{(m)}

with respect to

μ

. As seen from Eq. (28), the adjoint variables

λ_{s}^{(m)}

need to be computed in order to obtain the gradient. The adjoint variables are obtained from Eqs. (27) and (29) by starting with

m = M

. In this case,

Eq. (30)

{\tilde{W}}^{T} λ_{s}^{(M)} = {- (\frac{\partial f}{\partial {\tilde{Φ}}_{s}^{(M)}})}^{T} .

Then, for

m < M

the derivative

d f / d {\tilde{Φ}}_{s}^{(m)}

is obtained from Eq. (29) and used in Eq. (27) to get

λ_{s}^{(m)}

. From Eq. (27) it is seen that the matrix

{\tilde{W}}^{T}

needs to be repeatedly inverted in this process. The computation load of these inversions, which are always the same for a given algorithm iteration, can be reduced by first storing information on the matrix (e.g., some form of factorization). It should be pointed out that in solving the forward model, the inversion of

\tilde{W}

is needed as well. At that point, the information alluded to above regarding

\tilde{W}

can be stored.

As a final note, further applications of Eq. (29) are: 1. the construction of sensitivity maps (for each order $N$ ) to changes in the medium optical properties, and 2. the generation of single-perturbation error maps for custom data types and functional optimal selection.²²

4. Results of Numerical Experiments

4.1.

Reconstruction of Absorption Coefficient Maps

The numerical experiments in this subsection were performed for a circular geometry with a radius of 1.5 cm with homogeneous optical properties. This medium mimics a two dimensional (2D) slice of a small volume of biological tissue, see Fig. 2.

Fig. 2

Mesh of the circular geometry. Sources and detector (optodes) positions are represented by white spots (S) and dark squares (D), respectively.

The values of the optical coefficients of the homogeneous medium are $μ_{a} = 0.01 {cm}^{- 1}$ , $μ_{s} = 100 {cm}^{- 1}$ , and $g = 0.9$ . The refractive index of the medium is 1.4, resembling fat tissue, and the medium is considered as surrounded by air. These values are representative of small animal tissue samples in the NIR spectrum, e.g., abdomen, brain, and blood vessels in rats.²⁹^,⁵⁶ Next, Dirac delta functions in time collimated laser sources are placed at eight positions located at the boundary of the circular geometry, see Fig. 2. The sources are modeled as exponentially damped line sources distributed inside the medium and are fired sequentially. To use the FEM, the medium was divided into triangular elements (816 nodes and 1534 elements) using a coarse regular Delaunay triangulation as shown in Fig. 2. As a second step, an irregular-shaped inclusion is positioned inside the medium. It has the same properties as the background medium, except for the absorption coefficient $μ_{a}$ , for which the following increasing values were examined: 0.05, 0.1, and $1 {cm}^{- 1}$ (Fig. 3). These $μ_{a}$ values cover possible situations appearing in DOT and FDOT, including regions with high absorption, such as vascularized tissues and/or fluorescent inclusions.⁵⁶^,⁵⁷

Fig. 3

Circular geometry with an embedded absorptive inclusion (in white).

Using the FEM-FDM representation of the TD- $p {SP}_{N}$ equations, the time-dependent distribution of fluence in the circular geometry is calculated for the three aforementioned $μ_{a}$ values and for the orders $N = 1, 3, 5$ , and 7 of the ${SP}_{N}$ . A partially delimited darker region with depressed values of the fluence was noted at the location of the inclusion in all the cases. The darkening increased with the increase of the inclusion’s $μ_{a}$ . In addition, synthetic boundary data was generated by collecting the time-dependent values of the excitance. Measurements are gathered at eight different detector positions located at the boundary of the circular geometry (Fig. 2). Image reconstruction is then carried out with the DOT algorithm using the generated TD data and the same mesh used for the forward problem. At this point, it should be stressed that the inverse crime is strictly committed in the calculations (forward and inverse calculations employed the same mesh). This is done because the aimis to study the accuracy of the model (at each order $N$ ) in determining the corresponding set of optical parameters, without considering other sources of error, such as those introduced by using different mesh resolutions and/or noise in the measurements.To evaluate the image reconstructions, we introduce the relative error $ε_{r}$ (in percent) in the maximum of the reconstructed values at the region of interest (inclusion or background medium) $μ_{\max}^{re}$ , with respect to the original value $μ_{\max}^{o}$ as $ε_{r} = 100 \cdot | μ_{\max}^{re} - μ_{\max}^{o} | / μ_{\max}^{o}$ . The relative error $ε_{r}$ is also used in the multi-parameter reconstructions (next section) with the same purpose.

Figure 4 consecutively displays the solution to the inverse problem posed in Eq. (14) for the medium shown in Fig. 3. In all calculations, the absorption coefficient lower and upper bounds [see Eq. (14)] are taken to be 0.01 and $5 {cm}^{- 1}$ , respectively. For the three values of the inclusion’s $μ_{a}$ , 0.05, 0.1, and $1 {cm}^{- 1}$ (left to right), the inverse problem solution is plotted for the orders $N = 1$ , 3, and 5 (upper, middle, and lower rows, respectively). Order $N = 7$ was omitted from the exposed results because it did not introduce significant changes compared to order $N = 5$ . In all cases, the algorithm used the optical properties of the homogeneous medium as the initial guess; no a priori information was employed.

Fig. 4

Solution of the inverse problem for a homogeneous medium with an absorptive inclusion. The true values of the absorption coefficient for the inclusion are: 0.05 (1st column), 0.1 (2nd column), and $1 {cm}^{- 1}$ (3rd column). Absorption of the background is $0.01 {cm}^{- 1}$ . Solutions are plotted for the orders $N = 1$ , 3, and 5 (upper, middle, and lower rows, respectively). For comparison, the black ring located at the top of the figures represents the inclusion’s shape at a height corresponding to its absorption coefficient value. See text for precise numbers on the differences between reconstructed and true values.

Figure 4(a)–4(c) shows the reconstructed images using the ${SP}_{1}$ approximation can provide an accurate localization of the inclusion for the whole range of absorption coefficient values, and the background optical properties are accurately reconstructed; however, the absorption coefficient values of the heterogeneity are overestimated in all cases. The relative error $ε_{r}$ values range from 10% to 27% (left to right cases). This result is in accordance with the known fact that DE-based algorithms start failing in regions with high absorption values because the diffusion approximation on which the DE is based does not hold in such regions.³²

Figure 4(d)–4(f) shows the reconstructed images using the ${SP}_{3}$ approximation also accurately locate the position of the inclusion. The background optical properties are accurately recovered as well. The absorption coefficient values of the heterogeneity are underestimated in cases (d) and (f) and overestimated in case (e). The relative error $ε_{r}$ values vary from 5% (middle figure) to 9%, and13% (left and right figures).

Finally, Fig. 4(g)–4(i) shows the reconstructed images using the ${SP}_{5}$ approximation accurately locate the position of the inclusion and reconstruct the background optical properties. The reconstructed values of the inclusion’s absorption coefficient are overestimated in cases (h) and (i) and underestimated in case (g). The relative errors are 10%, 5%, and 19% (from left to right).

From the numerical experiments, it can be deduced that ${S P}_{N}$ orders higher than $N = 1$ (which corresponds to the DE) provide accurate localization of the absorptive heterogeneity and retrieval of the background optical properties. Moreover, they are able to better estimate high values of the absorption coefficient of the heterogeneity than the DE. Particularly, the reconstructed images using the ${SP}_{3}$ approximation provide the best estimates, with a reduction in the relative errors compared to the DE.

4.2.

Simultaneous Reconstruction of Absorption and Diffusion Maps

The simultaneous reconstruction of absorption and diffusion coefficient maps will now be discussed. The numerical experiments are performed with the same homogeneous circular geometry as previously (Fig. 2), and the optodes (sources and detectors) are in the same configuration as before. The optical coefficient values of the homogeneous medium are $μ_{a} = 0.01 {cm}^{- 1}$ , $μ_{s} = 80 {cm}^{- 1}$ , and $g = 0.9$ . These values correspond to a diffusion coefficient $D$ of 0.0416 cm. As before, the refractive index of the medium is taken as 1.4, and the medium is considered as surrounded by air. Besides pursuing the same goal as before, the inverse crime is again committed.

Next, two irregular inhomogeneities are embedded in the medium, see Fig. 5(a) and 5(c) (the inhomogeneities are shown in two different maps but they are really in the same medium; this is because the scales of absorption and scattering coefficients widely differ). The absorption coefficient of the first inclusion takes the increasing values of 0.05, 0.1, and $1 {cm}^{- 1}$ [Fig. 5(a) shows the case for $μ_{a} = 0.05 {cm}^{- 1}$ ] with remaining optical properties identical to those of the homogeneous background medium. These values correspond to diffusion coefficient values of 0.0414, 0.0412, and 0.037 cm, respectively.

The second inclusion [Fig. 5(c)] has the same properties as the background medium, except for the scattering coefficient $μ_{s}$ , which is assigned a value of $120 {cm}^{- 1}$ . This value corresponds to $D = 0.0277 cm$ . Note that the differences in the diffusion coefficient maps between the background and values at the inclusions are much higher for the scattering inclusion than the absorptive inclusion. As a consequence, for small absorption coefficient values, diffusion coefficient maps ( $D$ -maps) have to be carefully analyzed to appreciate the effect of absorptive inclusions. All optical coefficient values reported in this subsection are typical in small animal imaging (e.g., of commonly used in experimental sets).⁵⁸

As in the previous subsection, the time-dependent fluence distribution is calculated for all the physical situations and the orders $N = 1, 3, 5$ , and 7. Sources are modelled as previously. Differences in the fluence profiles partially delineate the zones where the inclusions are located. Excitance values are collected at eight detector positions located at the boundary of the circular geometry (Fig. 2). At this point, the absorption and the diffusion coefficient distributions are simultaneously reconstructed with the generated TD data.In the calculations, the absorption coefficient values [see Eq. (14)] are bound between 0.01 and $5 {cm}^{- 1}$ , which is the range of typical values encountered in practice. In a similar way, the diffusion coefficient is bound to lie between 0.0001 and 0.0416 cm (same value as the background medium), respectively. Reported error values for the optical parameters in this section are calculated as in the absorption-only reconstructions.

Fig. 5

(a) Absorption map depicting the absorptive inclusion ( $μ_{a} = 0.05 {cm}^{- 1}$ ). (b) Diffusion map (D-map) for the absorptive inclusion displayed in (a). (c) Scattering coefficient distribution maps for the scattering inclusion ( $μ_{s} = 120 {cm}^{- 1}$ ). (d) D-map for the scattering inclusion displayed in (c).

Figure 6 shows the reconstructed absorption coefficient maps ( $μ_{a}$ -maps). For the three values of $μ_{a}$ , 0.05, 0.1, and $1 {cm}^{- 1}$ (left to right columns of images), the inverse problem solution is plotted for the orders $N = 1, 3, 5$ , and 7 (first, second, third and fourth rows, respectively). In all cases, the algorithm uses the optical properties of the homogeneous medium as the initial guess; no a priori information is employed. In general, it can be seen that all orders ( $N = 1, 3, 5$ , and 7) provide an accurate localization of the inclusion for the whole range of $μ_{a}$ values presented. In addition, background optical properties are acceptably recovered, with an almost indistinguishable (less than 1%) relative error at few nodes (the reference used to calculate this error is the background absorption coefficient value).

Fig. 6

Solution of the inverse problem(absorption coefficient) for the multiparametric case. Values of the absorption coefficient are: 0.05, 0.1, and $1 {cm}^{- 1}$ (left to right columns of images). Images are plotted for the orders $N = 1$ , 3, 5, and 7 (first, second, third, and fourth rows, respectively). For comparison, the black ring located at the top of the figures represents the inclusion’s shape at a height corresponding to its absorption coefficient value. See text for precise numbers on the differences between reconstructed and true values.

For the absorptive inclusion with $μ_{a} = 0.05 {cm}^{- 1}$ (Fig. 6, first column of images), it is observed that the ${SP}_{1}$ and ${SP}_{7}$ approximations underestimate the value of $μ_{a}$ at the inclusion with errors of 19% and 12%, respectively. On the other hand, the ${SP}_{3}$ and ${SP}_{5}$ approximations overestimate the value of $μ_{a}$ at the inclusion with errors of 0.1% and 1%, respectively.

For $μ_{a} = 0.1 {cm}^{- 1}$ (Fig. 6, second column of images), the ${SP}_{1}$ approximation underestimates the value of $μ_{a}$ at the inclusion with a 16% error. In this physical situation, the orders $N = 3, 5$ , and 7 overestimate the value of $μ_{a}$ at the inclusion with errors of 8%, 12%, and 13%, respectively.

For $μ_{a} = 1 {cm}^{- 1}$ (Fig. 6, last column of images), only the ${SP}_{1}$ approximation overestimates the absorption coefficient value at the inclusion with an error of 8%. The ${SP}_{3}$ , ${SP}_{5}$ , and ${SP}_{7}$ approximations underestimate the absorption coefficient value at the inclusion, but with errors below 1%.

From the results, it may be inferred that in the case of multiparameter image reconstructions, orders higher than $N = 1$ provide an accurate localization of the absorptive heterogeneity and an acceptable quantitative retrieval of the background optical properties. In addition, they better quantify the inclusion’s $μ_{a}$ values than the DE. Reconstructed images using the ${SP}_{3}$ approximation give the best estimates, closely followed by the higher orders.

Figure 7 shows the reconstructed $D$ -maps (refer back to Fig. 6(b) and 6(d) for the true diffusion coefficient distribution). For the absorptive heterogeneity values of: 0.05, 0.1, and $1 {cm}^{1}$ (left to right), the inverse problem solution is plotted for the orders $N = 1, 3, 5$ , and 7 (first, second, third, and fourth rows, respectively). The diffusion coefficient reconstructions will now be described for each value of the absorption coefficient of the absorptive inclusion.

Fig. 7

Solution of the inverse problem(diffusion coefficient) for the multiparametric case. Images are plotted for the orders $N = 1$ , 3, 5, and 7 (first, second, third, and fourth rows, respectively) where the absorptive heterogeneity takes the following values: 0.05, 0.1, and $1 {cm}^{- 1}$ (left to right). These values are in turn associated with the following values of $D$ for the absorptive inclusion: 0.0414, 0.0412, and 0.0370 cm, respectively. As regards the scattering inclusion, the value of $D$ that should be recovered is 0.0277 cm. For comparison, a blue and a black ring are located in the figure to represent the absorption and the scattering inclusions’ shape at a height corresponding to its diffusion coefficient values. See text for precise numbers on the differences between reconstructed and true values.

For $μ_{a} = 0.05 {cm}^{- 1}$ (Fig. 7, first column of images), all ${SP}_{N}$ orders image reconstructions are able to accurately localize the scattering inclusion with a 5% to 6% error. The ${SP}_{3}$ gives better results, but only by a 1% better error than the others. The presence of the absorptive inclusion appears in the reconstructed image with a moderate intensity (in light gray in the image; only the peak of the distribution is clearly seen). Errors at this location are 9%, 7%, 9%, and 10% for the orders $N = 1, 3, 5$ , and 7, respectively. The background value of $D$ (0.0416 cm) is acceptably recovered by all the ${SP}_{N}$ orders. Some artifacts can be seen near the image boundary after an inspection of the image. Artifacts are more pronounced in the case of the ${SP}_{1}$ reconstructed image and in the high absorption case ( $1 {cm}^{- 1}$ ).

For $μ_{a} = 0.1 {cm}^{- 1}$ (Fig. 7, second column of images), all ${SP}_{N}$ image reconstructions are able to accurately localize the scattering inclusion. Errors for ${SP}_{1}$ , ${SP}_{3}$ , and ${SP}_{5}$ are the same: 6%. ${SP}_{7}$ shows a 20% error but with a better spatially delimited inclusion. As before, the presence of the absorptive inclusion appears in the reconstructed image with moderate intensity, although this time the visible spot is bigger due to a higher value of the absorptive inclusion. These results show the crosstalk effect between the absorption and diffusion coefficients and how it varies with the ${SP}_{N}$ order (the effect is strongest for $N = 7$ ). Errors at this location are 16%, 14%, 15%, and 3% for $N = 1, 3, 5$ , and 7, respectively. Small artifacts can be observed near the boundary, which again are more pronounced in the ${SP}_{1}$ reconstructed image.

For $μ_{a} = 1 {cm}^{- 1}$ (Fig. 7, last column of images), all ${SP}_{N}$ -based image reconstructions are able to accurately localize the scattering inclusion. However, the reconstructed images using orders $N = 3, 5$ , and 7 qualitatively provide a better shaped,recovered scattering inclusion in comparison with the original distribution. The errors in the inclusion’s scattering coefficient values are 6% for ${SP}_{3}$ and ${SP}_{5}$ and 15% for ${SP}_{7}$ . In the case of ${SP}_{1}$ reconstruction, the error is greater than 40%, showing its failure to accurately quantify the effect of the scattering inclusion in the presence of high absorption values. An examination of the presence of the absorptive inclusion in the $D$ -maps shows that the visible area has increased in size because of the large $μ_{a}$ value of the absorptive inclusion (cross-talk). The area is better delimited by the ${SP}_{3}$ and ${SP}_{7}$ reconstructions, followed by ${SP}_{5}$ reconstruction. Finally, ${SP}_{1}$ shows a highly irregular distribution at the absorptive inclusion. Errors at this location are difficult to quantify because of the irregular distribution of recovered optical coefficient values. The following estimates are obtained: 18%, 27%, and 8% for ${SP}_{3}$ , ${SP}_{5}$ , and ${SP}_{7}$ reconstructions, respectively. In the case of ${SP}_{1}$ reconstruction, the error is greater than 57%. As in the previous results, all images present small artifacts near the boundary, with more pronounced artifacts in the ${SP}_{1}$ reconstructed image.

From the previous results, it may be inferred that in the case of multiparameter image reconstructions, orders higher than $N = 1$ provide an accurate localization of the scattering heterogeneity and acceptable retrieval of the background optical properties. The ${SP}_{N}$ diffusion coefficient reconstructions provide images with better shaped inclusions for $N > 1$ when compared with the original distributions. Also, ${SP}_{N}$ ( $N > 1$ ) image reconstructions better quantify the diffusion coefficient values at the scattering inclusion than the DE reconstruction. In particular, the reconstructed images using the ${SP}_{3}$ approximation provide the best estimates for the diffusion coefficient values at the scattering inclusion. Furthermore, ${SP}_{N}$ ( $N > 1$ ) image reconstructions better quantify the absorptive heterogeneity contribution in the image than the DE. Particularly, ${SP}_{7}$ -based reconstructions give the best estimates for the diffusion coefficient values where the absorptive inclusion is located, but the size of the recovered inclusion is generally overestimated.

5. Comparison with Other Approaches

Theresults for multiparametric reconstructions obtained in the previous section will now be shortly compared with previously published works, Refs. 39, 42, 44, and 59. This is by no means an attempt to be exhaustive, but it will nevertheless give an idea about how our approach compares to others.

In Ref. 39, a DOT algorithm based on the FD ${SP}_{N}$ equations is presented. The forward model is discretized using the FEM (this builds on the ${SP}_{N}$ -based forward model presented in Ref. 31). The DOT algorithm uses an $L_{2}$ -norm Tikhonov minimization approach and a Levenberg-Marquardt procedure for updating the vector of optical properties. Multiparametric reconstructions (absorption and scattering coefficients) are performed with synthetic (noiseless) data. The inverse crime is strictly and purposely committed in the numerical experiments. The multi-parametric reconstructions appearing in this work can be appreciated in details in the presented figures. Comparing the results presented in Sec. 3.4 of Ref. 39 with those presented herein, it is found that inclusion optical properties are better estimated by the algorithm presented here. In addition, the algorithm fails to recover the scattering heterogeneity position when the order $N = 7$ is employed. The presence of cross-talk effects and image artifacts at the boundary is also mentioned by the authors. For $N = 7$ , image artifacts are located thorough the domain. Finally, high absorption regimes ( $\sim 1 {cm}^{- 1}$ ) are not studied in Ref. 39.

In Ref. 59, the time-dependent radiative transfer equation (RTE) is used in a model-based DOT algorithm. The RTE is spatially discretized on a structured grid and the angular dependence is approximated by a discrete ordinates procedure called the piecewise parabolic method (PPM), which ensures numerical stability.⁶⁰ Adjoint methods are employed in the calculation of the gradient of an objective square error function, similar to $f$ used here. The Polak–Ribière version of the nonlinear conjugate gradient method is used in the reconstructions, see details in the reference section of Ref. 59. 3D multiparameter reconstructions (absorption and scattering coefficients) using noiseless data are conducted. Although in two-dimensional (2D), the present algorithm provides better estimates of the inclusions optical properties for similar experimental cases. Hence, it may be expected that the present DOT algorithm can compete with that proposed in Ref. 59 in 3D cases. However, more conclusive experiments (same experimental conditions) would be needed, which is beyond the scope of the present work.

An analysis of the results presented in Refs. 42 and 44, where the RTE for FD problems is employed as the forward model and the inverse problem is posed as a PDE-constrained optimization, support the competitiveness of our DOT algorithm. In these works, the solution of the inverse problem is sought through an augmented Lagrangian⁴² and the reduced SQP⁴⁴ method. In Ref. 42, image artifacts and cross-talk effects which appear in the image reconstructions are difficult to evaluate for noise-free cases, mainly due to mesh effects (a different mesh is used for generating synthetic data and in the inverse problem solution). In Ref. 44, in the case of an absorptive inclusion and noise-free data, the background optical properties are not uniformly reconstructed, with large zones where the errors are around 20%. Finally, a considerable amount of image artifacts and cross-talk effects appears in the multiparametric reconstructions when 20 dB noise-added synthetic data is used. Hence, future works should be conducted to evaluate noise effects in our DOT algorithm and compared with these results in Ref. 44.

6. Conclusions

In this manuscript, a novel approach for solving inverse problems in time-domain diffuse optical tomography is presented. A finite element-finite difference discretized formulation of the time-dependent parabolic ${SP}_{N}$ equations is used as the forward model. In this way, complex geometries and heterogeneous distributions of the optical coefficients can be imaged. To achieve image reconstruction, an optimization method which includes bounds onthe optical coefficient values as additional linear constraints is implemented. The algorithm incorporates some necessary features for dealing with large-scale TD data. First, the implementation employs the nested analysis and design method to reduce the dimensionality of the problem. Therefore, even when TD data are dealt with, only the set of optical coefficients evaluated at the mesh nodes has to be recovered. Second, adjoint variables are introduced in the calculations. These reduce the computation time and provide an exact formula for the calculation of the gradient of the OF.

Several numerical experiments were performed for a small geometric medium with homogeneous background optical properties that mimic small animal imaging conditions. In a first round of experiments, an absorptive inclusion with increasing absorption coefficient values of 0.05, 0.1, and $1 {cm}^{- 1}$ was positioned into an otherwise homogeneous medium. These values were chosen to simulate practical situations encountered in DOT and FDOT small animal imaging. Synthetic TD data were generated in the simulations and used to solve the inverse problem for all the ${SP}_{N}$ orders studied in the literature. The inverse crime was committed in the numerical experiments, as the goal here was to evaluate the accuracy of each ${SP}_{N}$ order in determining the optical coefficient distributions. Noise-free TD data are employed in the calculations for the same reason. An analysis of the results showed that, in all the physical situations considered, the image reconstructions based on the ${SP}_{N}$ ( $N > 1$ ) provide an accurate localization of the absorptive heterogeneity and retrieval of the background optical properties. In addition, the ${SP}_{N}$ ( $N > 1$ ) better estimated the absorption coefficient values of the heterogeneity than the DE. In the simulations, the reconstructed images using the ${SP}_{3}$ approximation provide the best estimates of absorption coefficient distributions. The DE consistently provided the worst estimates with an error of near 20% for high absorption values at the inclusion.

In a second round of experiments, an absorptive inclusion with increasing absorption coefficient values of 0.05, 0.1, and $1 {cm}^{- 1}$ was placed in the homogeneous medium along with a second inclusion about half the size of the first, with optical properties identical to those of the homogeneous medium, except for the scattering coefficient to which a value of $120 {cm}^{- 1}$ was assigned. Using synthetic data collected at the boundary of the medium, a simultaneous reconstruction of the absorption and diffusion coefficient distributions was performed. The inverse crime was committed again in the numerical experiments. The absorption coefficient reconstructions showed that orders higher than $N = 1$ provide accurate localizations of the absorptive heterogeneity and an acceptable retrieval of the background optical properties. Moreover, they better quantify the absorption coefficient values of the inclusion than the DE. The ${SP}_{3}$ approximation stood outas the model providing the best results, followed by $N = 5$ and 7. The cross-talk effect appeared in the reconstructed images, reinforcing the absorption coefficient values of the inclusion. The cross-talk effect increases with increasing order $N$ . Similar results were gathered for the diffusion coefficient reconstructions. The image reconstructions based on the ${SP}_{N}$ ( $N > 1$ ) approximations provide accurate localizations of the scattering heterogeneity and acceptable retrieval of the background optical properties. In addition, they better quantify the values of the diffusion coefficient at the second inclusion (except for the ${SP}_{7}$ ), and the effect of the absorptive heterogeneity in the diffusion coefficient maps. Particularly, ${SP}_{3}$ approximation provides the best estimates for diffusion coefficient values at the scattering inclusion. The ${SP}_{7}$ approximation gives the best estimates for diffusion coefficient values at the absorptive inclusion when $μ_{a} = 0.1$ and $1 {cm}^{- 1}$ . The estimates of the optical properties of the absorptive inclusion by the DE reconstructions resulted in the worst values of all the orders, especially in the case of high absorption. In that case, the errors hamper any kind of estimation of the optical properties of scattering or absorptive inclusions. Cross-talk effects appeared in the reconstructed images in the form of small artifacts at the boundary, with these artifacts being more or less pronounced depending on the order. In the DE reconstructions, artifacts have a more significant contribution to the image.

The results obtained in this work indicate that our DOT algorithm can accurately substitute DE algorithms, especially in physical situations where the DE fails. The methods developed here can be used in further studies to 1. analyze the effect of other optical parameters in DOT reconstructions, such asthe anisotropy parameter, and perhaps eventually the refractive index, and 2. as a starting point for more elaborate optimization algorithms with other types offunctionals and additional constraints, and 3. it could eventually be reused with more accurate forward models.

Acknowledgments

J. Bouza-Domínguez acknowledges financial support from the FQRNT (Québec—Programme de bourses d’excellence pour étudiants étrangers—PBEEE). Y. Bérubé-Lauzière acknowledges financial support from an NSERC Discovery Grant (Canada) for the present work. Yves Bérubé-Lauzière is member of the FRQS-funded Centre de recherche clinique Étienne-Le Bel.

Appendices

Appendix A

In this Appendix, the structure of the matrices $C$ , $T$ , the matrix operator $D_{r}$ and the source vector $Q$ are provided up to the order $N = 7$ , which is the highest ${SP}_{N}$ order studied in the literature. In addition, the FEM matrices and vectors appearing in Eqs. (3) and (4) are given.

Matrix $C$ has the form (Matlab notation is used for specifying the columns)

Eq. (A1)

C (:, 1) = [\begin{matrix} μ_{0} (r) \\ - \frac{2}{3} μ_{0} (r) \\ \frac{8}{15} μ_{0} (r) \\ - \frac{16}{35} μ_{0} (r) \end{matrix}], C (:, 2) = [\begin{matrix} - \frac{2}{3} μ_{0} (r) \\ \frac{4}{9} μ_{0} (r) + \frac{5}{9} μ_{2} (r) \\ - \frac{16}{45} μ_{0} (r) - \frac{4}{9} μ_{2} (r) \\ \frac{32}{105} μ_{0} (r) + \frac{8}{21} μ_{2} (r) \end{matrix}], C (:, 3) = [\begin{matrix} \frac{8}{15} μ_{0} (r) \\ - \frac{16}{45} μ_{0} (r) - \frac{4}{9} μ_{2} (r) \\ \frac{64}{225} μ_{0} (r) + \frac{16}{45} μ_{2} (r) + \frac{9}{25} μ_{4} (r) \\ - \frac{128}{525} μ_{0} (r) - \frac{32}{105} μ_{2} (r) - \frac{54}{175} μ_{4} (r) \end{matrix}], C (:, 4) = [\begin{matrix} - \frac{16}{35} μ_{0} (r) \\ \frac{32}{105} μ_{0} (r) + \frac{8}{21} μ_{2} (r) \\ - \frac{128}{525} μ_{0} (r) - \frac{32}{105} μ_{2} (r) - \frac{54}{175} μ_{4} (r) \\ \frac{256}{1225} μ_{0} (r) + \frac{64}{245} μ_{2} (r) + \frac{324}{1225} μ_{4} (r) + \frac{13}{49} μ_{6} (r) \end{matrix}] .

The matrix

T

, and consequently

T^{- 1}

, do not depend on the medium optical coefficients and have the following structure

Eq. (A2)

T = [\begin{matrix} 1 & - \frac{2}{3} & \frac{8}{15} & - \frac{16}{35} \\ 0 & \frac{1}{3} & - \frac{4}{15} & \frac{8}{35} \\ 0 & 0 & \frac{1}{5} & - \frac{6}{35} \\ 0 & 0 & 0 & \frac{1}{7} \end{matrix}], T^{- 1} = [\begin{matrix} 1 & 2 & 0 & 0 \\ 0 & 3 & 4 & 0 \\ 0 & 0 & 5 & 6 \\ 0 & 0 & 0 & 7 \end{matrix}] .

The matrix operator

D_{r}

is a diagonal

4 \times 4

matrix operator whose diagonal elements are given by

Eq. (A3)

diag (D_{r}) = - [\begin{matrix} \nabla (\frac{1}{3 μ_{1}} \nabla), & \nabla (\frac{1}{7 μ_{3}} \nabla), & \nabla (\frac{1}{11 μ_{5}} \nabla), & \nabla (\frac{1}{15 μ_{7}} \nabla) \end{matrix}],

where the expression

diag ()

is used to list the diagonal elements. Finally, the source vector

Q (r, t)

contains the information about the source distribution

Q (r, t)

Eq. (A4)

Q (r, t) = Q (r, t) {[\begin{matrix} 1 & - \frac{2}{3} & \frac{8}{15} & - \frac{16}{35} \end{matrix}]}^{T} .

In the FEM, the domain of interest $V$ is partitioned into $l$ nonoverlapping elements $τ_{j}$ , $j = 1, \dots l$ joined at $d$ vertex nodes. The boundary of $V$ , denoted as $\partial V$ , will be assumed to have $d_{b}$ nodes. In the following, the index $l_{N} = (N + 1) / 2$ where $N$ is the order of the ${SP}_{N}$ approximation is used for convenience.

The square sparse matrix $\tilde{K}$ is a block diagonal matrix composed of elemental “stiffness” matrices ${{\tilde{K}}_{k}}_{k = 1, \dots, l_{N}}$ whose entries are given by

Eq. (A5)

{\tilde{K}}_{k} (i, j) = \int_{V} \frac{1}{(4 k - 1) μ_{2 k - 1}} \nabla u_{i} (r) \cdot \nabla u_{j} (r) d V, k = 1, \dots, l_{N}, i, j = 1, \dots, d,

where

μ_{k} = μ_{a} + μ_{s} (1 - g^{k})

are the transport coefficients, see Ref. 32. The square sparse matrices

\tilde{M}

,

\bar{Π}

, and

\tilde{T}

are composed of

l_{N} \times l_{N}

block “mass” matrices

{{\tilde{M}}_{k_{1}, k_{2}}}_{k_{1}, k_{2} = 1, \dots, l_{N}}

,

{{\tilde{Π}}_{k_{1}, k_{2}}}_{k_{1}, k_{2} = 1, \dots, l_{N}}

,

{{\tilde{T}}_{k_{1}, k_{2}}}_{k_{1}, k_{2} = 1, \dots, l_{N}}

whose entries are, respectively given by

Eq. (A6)

{\tilde{M}}_{k_{1}, k_{2}} (i, j) = \int_{V} C (k_{1}, k_{2}) u_{i} (r) u_{j} (r) d V, k_{1}, k_{2} = 1, \dots, l_{N}, i, j = 1, \dots, d,

Eq. (A7)

{\tilde{Π}}_{k_{1}, k_{2}} (i, j) = \int_{\partial V} \frac{Θ (k_{1}, k_{2})}{(4 k_{1} - 1) μ_{2 k_{1} - 1}} u_{i} (r) u_{j} (r) d σ, k_{1}, k_{2} = 1, \dots, l_{N}, i, j = 1, \dots, d,

Eq. (A8)

{\tilde{T}}_{k_{1}, k_{2}} (i, j) = \int_{V} T (k_{1}, k_{2}) u_{i} (r) u_{j} (r) d V, k_{1}, k_{2} = 1, \dots, l_{N}, i, j = 1, \dots, d .

In Eq. (A6) and (A8), $C (k_{1}, k_{2})$ and $T (k_{1}, k_{2})$ are the elements of the matrices $C$ and $T$ given in Eqs. (A1) and (A2), respectively. In Eq. (A7), the term $Θ = (k_{1}, k_{2})$ represents the elements of the matrix $Θ = {(B)}^{- 1} A$ which originates from the ${SP}_{N}$ boundary conditions, see Ref. 32. The matrices $A$ and $B$ have the following form

Eq. (A9)

A = [\begin{matrix} 1 / 2 + A_{1} & - 1 / 8 - C_{1} & 1 / 16 - E_{1} & - 5 / 128 - G_{1} \\ - 1 / 8 - C_{2} & 7 / 24 + A_{2} & - 41 / 384 - E_{2} & 1 / 16 - G_{2} \\ 1 / 16 - C_{3} & - 41 / 384 - E_{3} & 407 / 1920 + A_{3} & - 233 / 2560 - G_{3} \\ - 5 / 128 - C_{4} & 1 / 16 - E_{4} & - 233 / 2560 - G_{4} & 3023 / 17920 + A_{4} \end{matrix}],

Eq. (A10)

B = [\begin{matrix} (1 + B_{1}) / 3 μ_{1} & - D_{1} / μ_{3} & - F_{1} / μ_{5} & - H_{1} / μ_{7} \\ - D_{2} / 3 μ_{1} & (1 + B_{2}) / 7 μ_{3} & - F_{2} / μ_{5} & - H_{2} / μ_{7} \\ - D_{3} / 3 μ_{1} & - F_{3} / μ_{3} & (1 + B_{3}) / 11 μ_{5} & - H_{3} / μ_{7} \\ - D_{4} / 3 μ_{1} & - F_{4} / μ_{3} & - H_{4} / μ_{5} & (1 + B_{4}) / 15 μ_{7} \end{matrix}] .

The explicit expressions for the coefficients (

A_{1}, \dots, H_{1}

,

A_{4}, \dots, H_{4}

) can be found in Appendix A of Ref. 29. The vectors

{\tilde{F}}^{(m)}

and

{\tilde{Γ}}^{(m)}

are composed of

{{\tilde{F}}_{k}^{(m)}}_{k = 1, \dots, l_{N}}

and

{{\tilde{Γ}}_{k}^{(m)}}_{k = 1, \dots, l_{N}}

“elemental load vectors” evaluated at the time step

m

. Their entries are

Eq. (A11)

{\tilde{F}}_{k}^{(m)} (i) = \int_{V} Q^{(m)} (k) u_{i} (r) d V, k = 1, \dots, l_{N}, i = 1, \dots, d,

Eq. (A12)

{\tilde{Γ}}_{k}^{(m)} (i) = \int_{\partial V} \frac{G^{(m)} (k)}{(4 k - 1) μ_{2 k - 1}} u_{i} (r) d σ, k = 1, \dots, l_{N}, i = 1, \dots, d .

In Eq. (A11),

Q^{(m)} (k)

are the components of the source vector Eq. (A4) evaluated at time step

m

. In Eq. (A12),

G^{(m)} (k)

are the elements of the vector

G = B^{- 1} S

where the column vector

S

is related to the transmitted contribution of the exterior source (the exterior source power multiplied by the transmission coefficient)

B_{T} (r, \hat{s}, t)

as

Eq. (A13)

S = [\begin{matrix} \int_{\hat{s} \cdot \hat{n} > 0} B_{T} (r, \hat{s}, t) 2 | \hat{s} \cdot \hat{n} | d Ω \\ \int_{\hat{s} \cdot \hat{n} > 0} B_{T} (r, \hat{s}, t) [5 {| \hat{s} \cdot \hat{n} |}^{3} - 2 | \hat{s} \cdot \hat{n} |] d Ω \\ \int_{\hat{s} \cdot \hat{n} > 0} B_{T} (r, \hat{s}, t) [\frac{63}{4} {| \hat{s} \cdot \hat{n} |}^{5} - \frac{35}{2} {| \hat{s} \cdot \hat{n} |}^{3} + \frac{15}{4} | \hat{s} \cdot \hat{n} |] d Ω \\ \int_{\hat{s} \cdot \hat{n} > 0} B_{T} (r, \hat{s}, t) [\frac{429}{8} {| \hat{s} \cdot \hat{n} |}^{7} - \frac{693}{8} {| \hat{s} \cdot \hat{n} |}^{5} + \frac{315}{8} {| \hat{s} \cdot \hat{n} |}^{3} - \frac{35}{8} | \hat{s} \cdot \hat{n} |] d Ω \end{matrix}],

where

\hat{s}

and

\hat{n}

are unit vectors denoting direction and the normal to the boundary

\partial V

, respectively.

Appendix B

The vectors $j_{1}$ and $j_{2}$ have the following expressions

Eq. (B1)

j_{1} = [1 / 4 + J_{0}, (1 / 4 + J_{0}) (- 2 / 3) + (5 / 16 + J_{2}) (1 / 3), (1 / 4 + J_{0}) (8 / 15) + (5 / 16 + J_{2}) (- 4 / 15) + (- 3 / 32 + J_{4}) (1 / 5), \dots (1 / 4 + J_{0}) (- 16 / 35) + (5 / 16 + J_{2}) (8 / 35) + (- 3 / 32 + J_{4}) (- 6 / 35) + (13 / 256 + J_{6}) (1 / 7)],

Eq. (B2)

j_{2} = [- (\frac{0.5 + J_{1}}{3 μ_{1}}), (- \frac{J_{3}}{7 μ_{3}}), (- \frac{J_{5}}{11 μ_{5}}), (- \frac{J_{7}}{15 μ_{7}})],

where the coefficients

{J_{i}}_{i = 1, \dots, 4}

linearly depend on the angular moments of the angle-dependent Fresnel reflection coefficient

R_{F} (θ)

.⁶¹ The explicit expressions for the coefficients

J_{0}

,

J_{1, \dots,} J_{7}

can be found in Appendix A of Ref. 29.

References

1.

H. Dehghaniet al., “Numerical modelling and image reconstruction in diffuse optical tomography,” Phil. Trans. Royal Soc. A., 367 (1900), 3073 –3093 (2009). http://dx.doi.org/10.1098/rsta.2009.0090 PTRMAD 1364-503X Google Scholar

2.

S. R. ArridgeJ. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl., 25 (12), 123010 (2009). http://dx.doi.org/10.1088/0266-5611/25/12/123010 INPEEY 0266-5611 Google Scholar

3.

Translational Multimodality Optical Imaging, Artech House Inc, Boston, MA (2008). Google Scholar

4.

D. R. Leffet al., “Diffuse optical imaging of the healthy and diseased breast: a systematic review,” Breast Cancer Res. Treat., 108 (1), 9 –22 (2008). http://dx.doi.org/10.1007/s10549-007-9582-z BCTRD6 Google Scholar

5.

H. Zhaoet al., “Time-resolved diffuse optical tomography and its application to in vitro and in vivo imaging,” J. Biomed. Opt., 12 (6), 062107 (2007). http://dx.doi.org/10.1117/1.2815724 JBOPFO 1083-3668 Google Scholar

6.

J. C. Hebden, “Advances in optical imaging of the newborn infant brain,” Psychophysiology, 40 (4), 501 –510 (2003). http://dx.doi.org/10.1111/psyp.2003.40.issue-4 PSPHAF 0048-5772 Google Scholar

7.

A. P. GibsonJ. C. HebdenS. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol., 50 (4), R1 –R43 (2005). http://dx.doi.org/10.1088/0031-9155/50/4/R01 PHMBA7 0031-9155 Google Scholar

8.

R. ChoeA. CorluK. Lee, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys., 32 (4), 1128 –1139 (2005). http://dx.doi.org/10.1118/1.1869612 MPHYA6 0094-2405 Google Scholar

9.

Q. ZhangH. Jiang, “Three-dimensional diffuse optical tomography of simulated hand joints with a

64 \times 64

-channel photodiodes-based optical system,” J. Opt. A: Pure Appl. Opt., 7 (5), 224 –231 (2005). http://dx.doi.org/10.1088/1464-4258/7/5/003 JOAOF8 1464-4258 Google Scholar

10.

D. A. Boaset al., “Imaging the body with diffuse optical tomography,” IEEE Sign. Process. Mag., 18 (6), 57 –75 (2001). http://dx.doi.org/10.1109/79.962278 ISPRE6 1053-5888 Google Scholar

11.

A. Bluestoneet al., “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express, 9 (6), 272 –286 (2001). http://dx.doi.org/10.1364/OE.9.000272 OPEXFF 1094-4087 Google Scholar

12.

M. MarisE. GrattonJ. Maier, “Functional near-infrared imaging of deoxygenated hemoglobin during exercise of the finger extensor muscles using the frequency-domain technique,” Bioimaging, 2 (4), 174 –83 (1994). http://dx.doi.org/10.1002/1361-6374(199412)2:4<174::AID-BIO2>3.3.CO;2-H BOIMEL 0966-9051 Google Scholar

13.

D. K. Josephet al., “Diffuse optical tomography system to image brain activation with improved spatial resolution and validation with functional magnetic resonance imaging,” Appl. Opt., 45 (31), 8142 –8151 (2006). http://dx.doi.org/10.1364/AO.45.008142 APOPAI 0003-6935 Google Scholar

14.

J. Wanget al., “Spectral tomography with diffuse near-infrared light: inclusion of broadband frequency domain spectral data,” J. Biomed. Opt., 13 (4), 041305 (2008). http://dx.doi.org/10.1117/1.2952006 JBOPFO 1083-3668 Google Scholar

15.

H. Zhaoet al., “Time-resolved diffuse optical tomographic imaging for the provision of both anatomical and functional information about biological tissue,” Appl. Opt., 44 (10), 1905 –1916 (2005). http://dx.doi.org/10.1364/AO.44.001905 APOPAI 0003-6935 Google Scholar

16.

S. R. ArridgeW. R. B. Lionheart, “Nonuniqueness indiffusion-based optical tomography,” Opt. Lett., 23 (11), 882 –884 (1998). http://dx.doi.org/10.1364/OL.23.000882 OPLEDP 0146-9592 Google Scholar

17.

B. Harrach, “On uniqueness in diffuse optical tomography,” Inverse Probl., 25 (5), 055010 (2009). http://dx.doi.org/10.1088/0266-5611/25/5/055010 INPEEY 0266-5611 Google Scholar

18.

L. HervéA. KoenigJ.-M. Dinten, “Non-uniqueness in fluorescence-enhanced continuous wave diffuse optical tomography,” J. Opt., 13 (1), 015702 (2011). http://dx.doi.org/10.1088/2040-8978/13/1/015702 JOOPDB 0150-536X Google Scholar

19.

M. Schweigeret al., “Application of the finite element method for the forwardmodel in infrared absorption imaging,” Proc. SPIE, 1768 97 –108 (1992). http://dx.doi.org/10.1117/12.130893 PSISDG 0277-786X Google Scholar

20.

L. WangH. Wu, Biomedical Optics: Principles and Imaging, 1st ed.Wiley-Interscience, Hoboken, NJ (2007). Google Scholar

21.

F. Gaoet al., “The forward and inverse models in time-resolved optical tomography imaging and their finite-element method solutions,” Image Vision Comput., 16 (9–10), 703 –712 (1998). http://dx.doi.org/10.1016/S0262-8856(98)00078-X IVCODK 0262-8856 Google Scholar

22.

M. SchweigerS. R. Arridge, “Application of temporal filters to time-resolved data in optical tomography,” Phys. Med. Biol., 44 (7), 1699 –1717 (1999). http://dx.doi.org/10.1088/0031-9155/44/7/310 PHMBA7 0031-9155 Google Scholar

23.

F. E. W. Schmidtet al., “A 32-channel time-resolved instrument for medical optical tomography,” Rev. Sci. Instrum., 71 (1), 256 –265 (2000). http://dx.doi.org/10.1063/1.1150191 RSINAK 0034-6748 Google Scholar

24.

D. ColtonR. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, 133 2nd ed.Springer Verlag, New York (1998). Google Scholar

25.

F. GaoH. ZhaoY. Yamada, “Improvement of quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt., 41 (4), 778 –791 (2002). http://dx.doi.org/10.1364/AO.41.000778 APOPAI 0003-6935 Google Scholar

26.

H. Zhaoet al., “Time-resolved diffuse optical tomography and its application to in vitro and in vivo imaging,” J. Biomed. Opt., 12 (6), 062107 (2007). http://dx.doi.org/10.1117/1.2815724 JBOPFO 1083-3668 Google Scholar

27.

F. Gaoet al., “Simultaneous fluorescence yield and lifetime tomography from time-resolved transmittances of small-animal-sized phantom,” Appl. Opt., 49 (16), 3163 –3172 (2010). http://dx.doi.org/10.1364/AO.49.003163 APOPAI 0003-6935 Google Scholar

28.

A. H. HielscherR. E. AlcouffeR. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol., 43 (5), 1285 –1302 (1998). http://dx.doi.org/10.1088/0031-9155/43/5/017 PHMBA7 0031-9155 Google Scholar

29.

A. KloseE. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys., 220 (1), 441 –470 (2006). http://dx.doi.org/10.1016/j.jcp.2006.07.007 JCTPAH 0021-9991 Google Scholar

30.

Z. YuanH. Xin-HuaH. Jiang, “A higher order diffusion model for three-dimensional photon migration and image reconstruction in optical tomography,” Phys. Med. Biol., 54 (1), 65 –88 (2009). http://dx.doi.org/10.1088/0031-9155/54/1/005 PHMBA7 0031-9155 Google Scholar

31.

M. K. Chuet al., “Light transport in biological tissue using three-dimensional frequency-domain simplified spherical harmonics equations,” Phys. Med. Biol., 54 (8), 2493 –2509 (2009). http://dx.doi.org/10.1088/0031-9155/54/8/016 PHMBA7 0031-9155 Google Scholar

32.

J. Bouza DomínguezY. Bérubé-Lauzière, “Diffuse light propagation in biological media by a time-domain parabolic simplified spherical harmonics approximation with ray-divergence effects,” Appl. Opt., 49 (8), 1414 –1429 (2010). http://dx.doi.org/10.1364/AO.49.001414 APOPAI 0003-6935 Google Scholar

33.

A. LiemertA. Kienle, “Analytical solutions of the simplified spherical harmonics equations,” Opt. Lett., 35 (20), 3507 –3509 (2010). http://dx.doi.org/10.1364/OL.35.003507 OPLEDP 0146-9592 Google Scholar

34.

A. LiemertA. Kienle, “Comparison between radiative transfer theory and the simplified spherical harmonics theory for the semi-infinite geometry,” Opt. Lett., 36 (20), 4041 –4043 (2011). http://dx.doi.org/10.1364/OL.36.004041 OPLEDP 0146-9592 Google Scholar

35.

Y. Luet al., “A parallel adaptive finite element simplified spherical harmonics approximation solver for frequency domain fluorescence molecular imaging,” Phys. Med. Biol., 55 (16), 4625 –4645 (2010). http://dx.doi.org/10.1088/0031-9155/55/16/002 PHMBA7 0031-9155 Google Scholar

36.

J. Bouza DomínguezY. Bérubé-Lauzière, “Light propagation from fluorescent probes in biological tissues by coupled time-dependent parabolic simplified spherical harmonics equations,” Biomed. Opt. Express, 2 (4), 817 –837 (2011). http://dx.doi.org/10.1364/BOE.2.000817 BOEICL 2156-7085 Google Scholar

37.

L. D. MontejoH. K. KimA. H. Hielscher, “A finite volume algorithm for modeling light transport with the time-independent simplified spherical harmonics approximation to the equation of radiative transfer,” Proc. SPIE, 7896 78960J (2011). http://dx.doi.org/10.1117/12.875967 PSISDG 0277-786X Google Scholar

38.

Y. Luet al., “Spectrally resolved bioluminescence tomography with the third-order simplified spherical harmonics approximation,” Phys. Med. Biol., 54 (21), 6477 –6493 (2009). http://dx.doi.org/10.1088/0031-9155/54/21/003 PHMBA7 0031-9155 Google Scholar

39.

M. ChuH. Dehghani, “Image reconstruction in diffuse optical tomography based on simplified spherical harmonics approximation,” Opt. Express, 16 (12), 17780 –17791 (2009). http://dx.doi.org/10.1364/OE.17.024208 OPEXFF 1094-4087 Google Scholar

40.

J. Tianet al., “Evaluation of the simplified spherical harmonics approximation in bioluminescence tomography through heterogeneous mouse models,” Opt. Express, 18 (20), 20988 –21002 (2010). http://dx.doi.org/10.1364/OE.18.020988 OPEXFF 1094-4087 Google Scholar

41.

D. Hanet al., “Sparsity-promoting tomographic fluorescence imaging with simplified spherical harmonics approximation,” IEEE Tran. Biomed. Eng., 57 (10), 2564 –2567 (2010). http://dx.doi.org/10.1109/TBME.2010.2053538 IEBEAX 0018-9294 Google Scholar

42.

G. AbdoulaevK. RenA. H. Hielscher, “Optical tomography as a PDE-constrained optimization problem,” Inverse Probl., 21 (5), 1507 –153 (2005). http://dx.doi.org/10.1088/0266-5611/21/5/002 INPEEY 0266-5611 Google Scholar

43.

A. JoshiW. BangerthE. M. Sevick-Muraca, “Non-contact fluorescence optical tomography with scanning patterned illumination,” Opt. Express, 14 (14), 6516 –6534 (2006). http://dx.doi.org/10.1364/OE.14.006516 OPEXFF 1094-4087 Google Scholar

44.

H. K. KimA. H. Hielscher, “A PDE-constrained SQP algorithm for optical tomography based on the frequency-domain equation of radiative transfer,” Inverse Probl., 25 (1), 015010 (2009). http://dx.doi.org/10.1088/0266-5611/25/1/015010 INPEEY 0266-5611 Google Scholar

45.

H. K. Kimet al., “PDE-constrained multispectral imaging of tissue chromophores with the equation of radiative transfer,” Biomed. Opt. Express, 1 (3), 812 –824 (2010). http://dx.doi.org/10.1364/BOE.1.000812 BOEICL 2156-7085 Google Scholar

46.

H. K. KimJ. H. LeeA. H. Hielscher, “PDE-constrained fluorescence tomography with the frequency-domain equation of radiative transfer,” IEEE J. Sel. Top. Quant., 16 (4), 793 –803 (2010). http://dx.doi.org/10.1109/JSTQE.2009.2038112 IJSQEN 1077-260X Google Scholar

47.

S. B. Hazra, “PDE-constrained optimization methods,” Large-Scale PDE-Constrained Optimization in Applications, Lecture Notes in Applied and Computational Mechanics, 49 19 –25 Springer-Verlag, Berlin, Germany (2010). Google Scholar

48.

J. NocedalS. Wright, “Sequential quadratic programming,” Numerical Optimization, 529 –561 Springer, New York, NY (2006). Google Scholar

49.

A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, New York, NY (1978). Google Scholar

50.

A. V. PatilS. MukherjiU. B. Desai, “Optimal objective functional selection for image reconstruction indiffuse optical tomography,” in International Conference on Computing: Theory and Applications, 2007. ICCTA '07, 698 –704 (2007). Google Scholar

51.

N. Ducros, “Tomographie optique de fluorescence dans les milieux diffusants: apport de l’information temporelle,” Université Claude Bernard, (2009). Google Scholar

52.

F. Nouiziet al., “Improvement of absorption and scattering discrimination by selection of sensitive points on temporal profile in diffuse optical tomography,” Opt. Express, 19 (13), 12843 –12854 (2011). http://dx.doi.org/10.1364/OE.19.012843 OPEXFF 1094-4087 Google Scholar

53.

L. T. Biegleret al., “Numerical experience with a reduced Hessian method for large scale constrained optimization,” Comput. Optim. Appl., 15 (1), 45 –67 (2000). http://dx.doi.org/10.1023/A:1008723031056 CPPPEF 0926-6003 Google Scholar

54.

J. L. Huet al., “Sequential quadratic programming method for solution of electromagnetic inverse problems,” IEEE Trans. Antennas Propag., 53 (8), 2680 –2687 (2005). http://dx.doi.org/10.1109/TAP.2005.851871 IETPAK 0018-926X Google Scholar

55.

A. H. HielscherA. D. KloseK. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging, 18 (3), 262 –271 (1999). http://dx.doi.org/10.1109/42.764902 ITMID4 0278-0062 Google Scholar

56.

D. C. ComsaT. J. FarrellM. S. Patterson, “Quantitative fluorescence imaging of point-like sources in small animals,” Phys. Med. Biol., 53 (20), 5797 –5814 (2008). http://dx.doi.org/10.1088/0031-9155/53/20/016 PHMBA7 0031-9155 Google Scholar

57.

J. MobleyT. Vo-Dinh, “Optical properties of tissue,” Biomedical Photonics Handbook, 20 –95 CRC Press, Boca Raton, FL (2003). Google Scholar

58.

G. AlexandrakisF. R. RannouA. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol., 50 (17), 4225 –4241 (2005). http://dx.doi.org/10.1088/0031-9155/50/17/021 PHMBA7 0031-9155 Google Scholar

59.

J. BoulangerA. Charette, “Numerical developments for short-pulsednear-infrared spectroscopy. Part II: inverse treatment,” JQSRT, 91 (3), 297 –318 (2005). http://dx.doi.org/10.1016/j.jqsrt.2004.05.062 JQSRAE 0022-4073 Google Scholar

60.

J. BoulangerA. Charette, “Numerical developments for short-pulsed near-infrared spectroscopy. Part I: direct treatment,” JQSRT, 91 (2), 189 –209 (2005). http://dx.doi.org/10.1016/j.jqsrt.2004.05.057 JQSRAE 0022-4073 Google Scholar

61.

M. BornE. Wolf, Principles of Optics, 7th ed.Cambridge University Press, Cambridge, UK (2003). Google Scholar

Citation Download Citation

Jorge Bouza Domìnguez and Yves Bérubé-Lauzière "Diffuse optical tomographic imaging of biological media by time-dependent parabolic SPN equations: a two-dimensional study," Journal of Biomedical Optics 17(8), 086012 (29 August 2012). https://doi.org/10.1117/1.JBO.17.8.086012

Published: 29 August 2012

Access the abstract

JOURNAL ARTICLE
8 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY