Nonlinear approach to difference imaging in diffuse optical tomography

Meghdoot Mozumder; Tanja Tarvainen; Aku Seppänen; Ilkka Nissilä; Simon R. Arridge; Ville Kolehmainen

doi:10.1117/1.JBO.20.10.105001

6 October 2015 Nonlinear approach to difference imaging in diffuse optical tomography

Meghdoot Mozumder, Tanja Tarvainen, Aku Seppänen, Ilkka Nissilä, Simon R. Arridge, Ville Kolehmainen

Author Affiliations +

Journal of Biomedical Optics, Vol. 20, Issue 10, 105001 (October 2015). https://doi.org/10.1117/1.JBO.20.10.105001

Abstract

Difference imaging aims at recovery of the change in the optical properties of a body based on measurements before and after the change. Conventionally, the image reconstruction is based on using difference of the measurements and a linear approximation of the observation model. One of the main benefits of the linearized difference reconstruction is that the approach has a good tolerance to modeling errors, which cancel out partially in the subtraction of the measurements. However, a drawback of the approach is that the difference images are usually only qualitative in nature and their spatial resolution can be weak because they rely on the global linearization of the nonlinear observation model. To overcome the limitations of the linear approach, we investigate a nonlinear approach for difference imaging where the images of the optical parameters before and after the change are reconstructed simultaneously based on the two datasets. We tested the feasibility of the method with simulations and experimental data from a phantom and studied how the approach tolerates modeling errors like domain truncation, optode coupling errors, and domain shape errors.

1. Introduction

Diffuse optical tomography (DOT) is a noninvasive technique for imaging biological tissues with applications in imaging of human brain¹^–³ and breast⁴^,⁵ and small animals.⁶^,⁷ In brain imaging, DOT has been used for functional brain activation studies,¹^,³^,⁸^–¹⁰ imaging hemorrhages at birth,² and stroke.⁷^,¹¹ DOT is portable, which makes it a potential imaging method for bed-side monitoring of newborn infants or adults in intensive care.⁸ Recently, it has also been suggested that using high-density imaging arrays DOT could achieve spatial resolution comparable to functional magnetic resonance imaging.¹²^,¹³ DOT has also applications in cancer tumor diagnosis in humans and small animals.⁴^–⁶

Absolute imaging in DOT uses a single set of measurements to reconstruct spatially distributed absorption and scattering coefficients. However, absolute imaging is very sensitive to modeling errors, which can be caused, e.g., by inaccurately known object shape or unknown optode coupling coefficients or optode positions. A variety of techniques have been developed that could provide tolerance toward such modeling errors. Some of the techniques rely on explicit calculation of coupling coefficients from experiments¹⁴^,¹⁵ or using computational methods.¹⁶^,¹⁷ For estimating the domain shape, several registration methods are available that typically use a finite set of measured surface points to fit a generic head model, in the form of an anatomical atlas.¹⁸ However, interpolating the object shape using a few measured points does not guarantee obtaining the exact surface of the patient; hence, the process might still retain modeling errors. Also, techniques using data from other imaging modalities such as computed tomography (CT)¹⁹ or magnetic resonance imaging (MRI)²⁰^,²¹ to obtain the domain shape and optode positions have been proposed. Such data might not always be available. The Bayesian approximation error approach²² is an alternative computational technique where the statistics of such model-based errors are precomputed using prior probability distributions of the unknowns and the uncertain nuisance parameters. These statistics are then used in the image reconstruction process to compensate for the modeling errors.²³^–²⁸

Nevertheless, the most popular technique for in vivo DOT imaging has been difference imaging where the objective is to reconstruct change in the optical properties using measurements before and after the change. Conventionally, the image reconstruction is carried out using the difference of the measurements and a linearized approximation of the observation model.³^,⁸^–¹⁰^,¹⁴^,²⁹^–³² In the case of imaging brain activation, the reference measurement from the state before the change is typically obtained at the “rest stage” of a brain.⁸^,⁹ In some cases, the reference measurement can also be obtained from a homogeneous tissue mimicking phantom.³² In Refs. 29 30.–31, 33, the difference imaging problem is stated as reconstructing moment-by-moment relative differences of the parameters using mean of the time series of measured data as the reference data. One of the main benefits of the linearized difference imaging approach is that it has a good tolerance to (invariant) modeling errors, such as inaccurately known source and detector locations and coupling as well as inaccurately known body shape. When reconstructing images using differences in data resulting from a change in the optical properties, several modeling errors, which are invariant between the measurements, cancel out (partially) in the subtraction of the measurement before the change from the measurement after the change.¹⁴ The performance of linear difference imaging in the presence of mismodeled background was studied in Refs. 25, 31 and in imaging objects of different sizes and target optical properties in Ref. 29. The method was extended for extracting dynamic information from time series data in Refs. 33 and 34. A nonlinear (postprocessing) update method based on spatial deconvolution of the difference images was implemented in Refs. 34 35.36.37.38.–39 to improve on the linearized solutions. A drawback of the linear reconstruction approach is that the difference images are usually only qualitative in nature and their spatial resolution can be weak because they rely on the global linearization of the nonlinear observation model. The performance of the linear reconstruction also depends on the linearization point, which ideally should be equivalent to the initial state, which in practice is always unknown. According to simulations in Ref. 40, inaccurately known background optical properties together with linear difference imaging can lead to inaccurate contrast in the reconstructed images, and in some cases it is possible that the method even fails to detect and localize the changes.

To overcome the limitations of the linear reconstruction approach, we apply a new nonlinear approach for difference imaging in DOT that has been initially developed for estimating conductivity distribution in electrical impedance tomography in Refs. 41 and 42. In this approach, the optical parameters after the change are parameterized as a linear combination of the initial state and the change. The approach is based on the regularized nonlinear least-squares approach. Instead of using the difference of the data before and after the change, the measurement before and after the change is concatenated into a single measurement vector, and the objective is to estimate the unknown initial state and the change based on the data. This approach allows naturally for modeling independently the spatial characteristics of the background optical parameters and the change of the optical parameters by separate regularization functionals. The approach also allows the restriction of the optical parameter changes to a region of interest (ROI) inside the domain in cases where the change is known to occur in a certain subvolume of the body. We test the feasibility of the method with two-dimensional (2-D) simulations using frequency domain data from a simplified head geometry and experimental frequency domain data from a cylindrical phantom using three-dimensional (3-D) models.

Since good tolerance for modeling errors such as domain truncation, inaccurately known optode coupling coefficients, and inaccurately known domain shape are some of the main benefits of linear difference imaging, we also study how the new nonlinear approach tolerates the same modeling errors.

The remainder of the paper is organized as follows. In Sec. 2, a brief review of the light transport modeling in DOT is given. Next, the absolute and difference imaging approaches in DOT and the reconstruction algorithms are explained. In Sec. 3, we describe the methods used in data simulation and constructing the prior models. Then, we present the results in a 2-D geometry with and without modeling errors. Next, we describe our experimental setup and present the reconstruction results with the experimental data. Finally, the conclusions are given in Sec. 4.

2. Forward Model of Diffuse Optical Tomography

Let $Ω \subset R^{n}$ , $n = 2, 3$ , denote the object domain and $\partial Ω$ the domain boundary. In a diffusive medium, the commonly used light transport model for DOT is the diffusion approximation to the radiative transport equation. ⁴³ In this paper, the frequency domain version of the diffusion approximation is used,

Eq. (1)

(- \nabla \cdot κ (r) \nabla + μ_{a} (r) + \frac{j ω}{c}) Φ (r) = 0, r \in Ω,

Eq. (2)

Φ (r) + \frac{1}{2 γ} κ (r) α \frac{\partial Φ (r)}{\partial \hat{k}} = {\begin{cases} \frac{Q_{p}}{γ} & r \in s_{p} \\ 0 & r \in \partial Ω \ s_{p} \end{cases},

where

Φ (r) = Φ

is the fluence,

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

is the absorption coefficient, and

κ (r) = κ

is the diffusion coefficient. The diffusion coefficient

κ

is given by

κ (r) = 1 / {n [μ_{a} (r) + {μ_{s}}^{'} (r)]}

, where

μ_{s}^{'} (r) = μ_{s}^{'}

is the (reduced) scattering coefficient. Furthermore,

j

is the imaginary unit,

ω

is the angular modulation frequency of the input signal, and

c

is the speed of light in the medium. The parameter

Q_{p}

is the strength of the light source at locations

s_{p}

,

p = 1, \dots, N_{s} \subset \partial Ω

, operating at angular frequency

ω

. The parameter

γ

is a dimension-dependent constant (

γ = 1 / π

when

Ω \subset R^{2}

,

γ = 1 / 4

when

Ω \subset R^{3}

) and

α

is a parameter governing the internal reflection at the boundary

\partial Ω

. The measurable quantity exitance

Γ (r)

by detector

q

under illumination from source

p

is given by

Eq. (3)

Γ_{p q} (r) = - κ \frac{\partial Φ_{p} (r)}{\partial \hat{k}} = \frac{2 γ}{α} Φ_{p} (r) r \in d_{q},

where

\hat{k}

is the outward normal to the boundary at point

r

and

d_{q}

,

q = 1, \dots, N_{d} \subset \partial Ω

are the detector locations.

The numerical approximation of the forward model (1)–(3) is based on a finite-element (FE) approximation. In the FE approximation, the domain $Ω$ is divided into $N_{e}$ nonoverlapping elements joined at $N_{n}$ vertex nodes. The photon density in the finite dimensional basis is given by

Eq. (4)

Φ^{h} = \sum_{k = 1}^{N_{n}} φ_{k} ψ_{k} (r),

where

ψ_{k}

is the nodal basis functions of the FE mesh and

φ_{k}

is the photon densities in the nodes of the FE mesh. Furthermore, we write finite dimensional (piecewise constant) approximations for

μ_{a} (r)

and

μ_{s}^{'} (r)

Eq. (5)

μ_{a} (r) = \sum_{l = 1}^{N} μ_{a, l} χ_{l} (r),

Eq. (6)

μ_{s}^{'} (r) = \sum_{l = 1}^{N} μ_{s, l} χ_{l} (r),

where

χ_{l}

denotes the characteristic functions of disjoint image pixels.

The measurement data for frequency domain DOT typically contains the measured log amplitude and phase for all source-detector pairs

Eq. (7)

y = (\begin{matrix} Re \log (Γ) \\ Im \log (Γ) \end{matrix}) \in R^{2 N_{s} N_{d}},

where

y

is the data vector. The FE-based solution of Eqs. (1)–(3) is denoted by

A (x)

, where

x = {(μ_{a}, μ_{s}^{'})}^{T} \in R^{2 N}

is the discretized optical coefficients. The observation model is written as

Eq. (8)

y = A (x) + e,

where

e \in R^{2 N_{s} N_{d}}

models the random noise in measurements.

2.1.

Absolute Imaging

In absolute imaging, the optical coefficients $x$ are reconstructed using a single set of measurements $y$ during which the target is assumed to be nonvarying.

Assuming that the additive measurement noise is independent of the unknowns and distributed as zero-mean Gaussian $e \sim N (0, Γ_{e})$ ,⁴⁴ the estimate amounts to the minimization problem

Eq. (9)

\hat{x} = \arg \min_{x > 0} {{∥ L_{e} [y - A (x)] ∥}^{2} + f (x)},

where the

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

is the Cholesky factor and

f (x)

is the regularization functional, which should be constructed based on the prior information on the unknowns.

Usual choices for the regularization functional include standard Tikhonov regularization $f (x) = α {∥ x ∥}^{2}$ , smoothness regularization $f (x) = α {∥ L_{x} x ∥}^{2}$ , where $L_{x}$ is a (possibly spatially and directionally weighted) differential operator,³^,⁴⁵^,⁴⁶ total variation (TV)⁴⁷^–⁵⁰ regularization $f (x) = α {∥ \nabla x ∥}_{1}$ , and so on.

We note that the estimate (9) can be interpreted in the Bayesian inversion framework as the maximum a posteriori estimate from a posterior density model, which is based on the observation model (8) and a prior model for the unknowns.²²^,²³^,⁵¹

2.2.

Difference Imaging

Consider two DOT measurement realizations $y_{1}$ and $y_{2}$ obtained from the body at times $t_{1}$ and $t_{2}$ with optical coefficients $x_{1}$ and $x_{2}$ , respectively. The observation models corresponding to the two DOT measurement realizations can be written as in Eq. (8)

Eq. (10)

y_{1} = A (x_{1}) + e_{1},

Eq. (11)

y_{2} = A (x_{2}) + e_{2},

where

e_{i} \sim N (0, Γ_{e_{i}})

,

i = 1,2

. The aim in difference imaging is to reconstruct the change in optical parameters

δ x = x_{2} - x_{1}

based on the measurements

y_{1}

and

y_{2}

.

2.2.1.

Conventional linear approach to difference imaging

Conventionally, the image reconstruction in difference imaging is carried out as follows. Equations (10) and (11) are approximated by the first-order Taylor approximations as

Eq. (12)

y_{i} \approx A (x_{0}) + J (x_{i} - x_{0}) + e_{i}, i = 1,2,

where

x_{0}

is the linearization point and

J = \partial A / \partial x (x_{0})

is the Jacobian matrix evaluated at

x_{0}

. Using the linearization and subtracting

y_{1}

from

y_{2}

gives the linear observation model

Eq. (13)

δ y \approx J δ x + δ e,

where

δ x = x_{2} - x_{1}

and

δ e = e_{2} - e_{1}

.

Given the model (13), the change in optical coefficients $δ x$ can be estimated as

Eq. (14)

\hat{δ x} = \arg \min_{δ x} {{‖ L_{δ e} (δ y - J δ x) ‖}^{2} + f_{δ x} (δ x)},

where

f_{δ x} (δ x)

is the regularization functional. The weighting matrix

L_{δ e}

is defined as

L_{δ e}^{T} L_{δ e} = Γ_{δ e}^{- 1}

, where

Γ_{δ e} = Γ_{e_{1}} + Γ_{e_{2}}

.

The regularization functional $f_{δ x} (δ x)$ is often chosen to be of the quadratic form $f_{δ x} (δ x) = {‖ L_{δ x} δ x ‖}^{2}$ , where $L_{δ x}$ is the regularization matrix. In such a case, the problem (14) is linear and has a closed form solution¹³^,⁵²

Eq. (15)

\hat{δ x} = {(J^{T} Γ_{δ e}^{- 1} J + Γ_{δ x}^{- 1})}^{- 1} (J^{T} Γ_{δ e}^{- 1} δ y) .

In this paper, we refer to Eq. (15) as the conventional difference imaging estimate.

The main benefit of the difference imaging is that at least part of the modeling errors cancel out when considering the difference data $δ y$ . Hence, the estimates are often to some extent tolerant of modeling errors. A drawback of the approach is that the difference images are usually only qualitative in nature and their spatial resolution can be weak because they rely on the global linearization of the nonlinear observation model (8). Moreover, the estimates depend on the selection of the linearization point $x_{0}$ . Typically, $x_{0}$ is selected as a homogeneous (spatially constant) estimate of the initial state $x_{1}$ . This choice can lead to errors in the reconstructions if the initial optical coefficients are not accurately known.²⁵^,³¹

2.2.2.

Nonlinear approach to difference imaging

In this section, we formulate the reconstruction of the change of optical coefficients in the nonlinear regularized least squares framework. However, instead of reconstructing both states $x_{1}$ and $x_{2}$ using Eq. (9) separately for datasets $y_{i}$ , $i = 1,2$ and then subtracting $δ x = x_{2} - x_{1}$ , we use an approach where $δ x$ is reconstructed together with $x_{1}$ by using both datasets $y_{1}$ and $y_{2}$ simultaneously.⁴¹^,⁴² This approach allows us to model prior information in cases where, e.g., the spatial characteristics of the initial state $x_{1}$ and change $δ x$ are different (e.g., smooth $x_{1}$ and sparse $δ x$ ). This approach also allows in a straightforward way the restriction of the change $δ x$ into a ROI.

Let us assume that the change in optical coefficients $δ x = x_{2} - x_{1}$ is known to occur in $Ω_{ROI} \subseteq Ω$ , and denote the change of optical coefficients within $Ω_{ROI}$ by $δ x_{ROI}$ . Then, $δ x = M δ x_{ROI}$ , where $M$ is an extension mapping $M : Ω_{ROI} \to Ω$ such that

Eq. (16)

M δ x_{ROI} (r) = {\begin{cases} δ x_{ROI} (r), & r \in Ω_{ROI} \\ 0 & r \in Ω \ Ω_{ROI} \end{cases} .

Obviously, if no ROI constraint is used, we set $Ω_{ROI} = Ω$ . The optical coefficients after the change $x_{2}$ can now be represented as a linear combination of the initial state and the change as

Eq. (17)

x_{2} = x_{1} + M δ x_{ROI} .

Inserting Eq. (17) into Eq. (11) and concatenating the measurement vectors $y_{1}$ and $y_{2}$ and the corresponding models in Eqs. (10) and (11) into block vectors, leads to an observation model⁴¹

Eq. (18)

\binom{\underset{⏟}{[\begin{matrix} y_{1} \\ y_{2} \end{matrix}]}}{\tilde{y}} = \binom{\underset{⏟}{[\begin{matrix} A (x_{1}) \\ A (x_{1} + M δ x_{ROI}) \end{matrix}]}}{\tilde{A} (\tilde{x})} + \binom{\underset{⏟}{[\begin{matrix} e_{1} \\ e_{2} \end{matrix}]}}{\tilde{e}},

Eq. (19)

\tilde{y} = \tilde{A} (\tilde{x}) + \tilde{e},

where

Eq. (20)

\tilde{x} = [\begin{matrix} x_{1} \\ δ x_{ROI} \end{matrix}] .

Given the model in Eq. (19), the initial state $x_{1}$ and the change $δ x$ can be simultaneously estimated as

Eq. (21)

\hat{\tilde{x}} = \arg \min_{x_{1} > 0} {{∥ L_{\tilde{e}} [\tilde{y} - \tilde{A} (\tilde{x})] ∥}^{2} + f (\tilde{x})} .

Here, $L_{\tilde{e}} \in R^{4 N_{s} N_{d} \times 4 N_{s} N_{d}}$ is the Cholesky factor such that $L_{\tilde{e}}^{T} L_{\tilde{e}} = Γ_{\tilde{e}}^{- 1}$ , where

Γ_{\tilde{e}} = [\begin{matrix} Γ_{e_{1}} & 0_{2 N_{s} N_{d} \times 2 N_{s} N_{d}} \\ 0_{2 N_{s} N_{d} \times 2 N_{s} N_{d}} & Γ_{e_{2}} \end{matrix}],

and

0_{2 N_{s} N_{d} \times 2 N_{s} N_{d}} \in R^{2 N_{s} N_{d} \times 2 N_{s} N_{d}}

is an all-zero matrix.

f (\tilde{x})

is the joint regularization functional of

\tilde{x} = {(x_{1}, δ x)}^{T}

, which allows for separate models for

x_{1}

and

δ x

as

Eq. (22)

f (\tilde{x}) = f_{1} (x_{1}) + f_{2} (δ x) .

The estimate in Eq. (21) can be computed iteratively using for example a Gauss–Newton (GN) algorithm⁵³ as

Eq. (23)

{\hat{\tilde{x}}}_{i + 1} = {\hat{\tilde{x}}}_{i} + {\hat{s}}_{i} {[{\tilde{J}}_{i}^{T} Γ_{\tilde{e}}^{- 1} {\tilde{J}}_{i} + ℋ ({\hat{\tilde{x}}}_{i})]}^{- 1} [{\tilde{J}}_{i}^{T} Γ_{\tilde{e}}^{- 1} {\tilde{y} - \tilde{A} ({\hat{\tilde{x}}}_{i})] - G ({\hat{\tilde{x}}}_{i})},

where

{\hat{s}}_{i}

is the step length, and

G ({\hat{\tilde{x}}}_{i})

and

ℋ ({\hat{\tilde{x}}}_{i})

are the gradient and Hessian of the regularization functional

f (\tilde{x})

, respectively. The Jacobian matrix

{\tilde{J}}_{i} = \partial \tilde{A} / \partial \tilde{x}

is of the form

{\tilde{J}}_{i} (\hat{\tilde{x}}) = [\begin{matrix} J ({\hat{x}}_{1, i}) & 0_{2 N_{s} N_{d} \times N_{ROI}} \\ J ({\hat{x}}_{1, i} + M \hat{δ x_{ROI, i}}) & J ({\hat{x}}_{1, i} + M \hat{δ x_{ROI, i}}) M \end{matrix}],

where

J (x)

is the Jacobian matrix for the forward mapping

A (x)

,

0_{2 N_{s} N_{d} \times N_{ROI}} \in R^{2 N_{s} N_{d} \times N_{ROI}}

is an all-zero matrix, and

N_{ROI}

is the dimension of the vector

δ x_{ROI}

.

3. Results

The feasibility of the method was tested with 2-D simulations and with 3-D experimental measurement data.

3.1.

Estimates

The following reconstruction approaches were considered:

Nonlinear: nonlinear reconstruction of $\tilde{x} = {(x_{1}, δ x_{ROI})}^{T}$ with the proposed (nonlinear) difference imaging method

Eq. (24)

\hat{\tilde{x}} = \arg \min_{x_{1} > 0} {{∥ L_{\tilde{e}} [\tilde{y} - \tilde{A} (\tilde{x})] ∥}^{2} + f_{1} (x_{1}) + f_{2} (δ x_{ROI})} .

Linear: conventional (linear) difference reconstruction by solving

Eq. (25)

\hat{δ x} = \arg \min_{δ x} {{∥ L_{δ e} (δ y - J δ x) ∥}^{2} + f (δ x)} .

3.2.

Regularization Functionals

For modeling $x_{1}$ in the nonlinear difference imaging approach Eq. (24), we used a quadratic Gaussian smoothness regularization of the form

Eq. (26)

f_{1} (x_{1}) = {∥ L_{x_{1}} (x_{1} - x_{1, *}) ∥}^{2} .

In the construction of the Gaussian smoothness regularization functional Eq. (26), the absorption and scatter coefficients $μ_{a, 1}$ and $μ_{s, 1}^{'}$ were modeled as mutually independent Gaussian random fields. In the construction of prior covariances, the random field $x_{1}$ was considered in the form

x_{1} = x_{1, in} + x_{1, bg},

where

x_{1, in}

is a spatially inhomogeneous parameter with zero mean,

x_{1, in} \sim N (0, Γ_{in, x_{1}})

, and

x_{1, bg}

is a spatially constant (background) parameter with nonzero mean. For the latter, we can write

x_{bg} = I

, where

I

is a vector of ones and

q

is a scalar random variable with Gaussian distribution

q

\sim N (c, σ_{bg, x_{1}}^{2})

. In the construction of

Γ_{in, f}

, the approximate correlation length was adjusted to match the expected size of the inhomogeneities and the marginal variances of

x_{in}

were tuned based on the expected variation of the optical properties in the initial state. See Refs. 22, 23, and 26 for further details. Modeling

x_{1, in}

and

x_{1, bg}

as mutually independent, we obtain

x_{1, *} = c I, Γ_{x_{1}} = Γ_{in, x_{1}} + σ_{bg, x_{1}}^{2} {I I}^{T},

where the Cholesky factor

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

in Eq. (26). This particular Gaussian smoothness regularization has been previously used in absolute imaging in Refs. 23 24.–25, 27, 28, 48 and in difference imaging in Ref. 54. In this paper, the correlation length was chosen as 8 mm. The standard deviations of the background and inhomogeneities are given in Table 1.

Table 1

Parameters of Gaussian smoothness regularization. The first two rows show standard deviations of the background σbg,δx and the inhomogeneities σin,δx chosen for modeling δx. The next two rows show σbg,x1 and σin,x1 chosen for modeling the initial state x1. Here, x1,* is selected as a homogeneous (spatially constant) estimate of the initial state x1.

Regularization parameters	(2-D and 3-D)
$σ_{bg, δ x}$	$0.1 \times x_{1, *}$
$σ_{in, δ x}$	$0.2 \times x_{1, *}$
$σ_{bg, x_{1}}$	$0.4 \times x_{1, *}$
$σ_{in, x_{1}}$	$0.8 \times x_{1, *}$

For modeling $δ x_{ROI}$ (in all the following test cases except for case 2 in Sec. 3.3.2), we used a sparsity promoting TV functional $f_{2} (δ x_{ROI}) = α TV (δ x_{ROI})$ , where

Eq. (27)

TV (δ x_{ROI}) = \sum_{k = 1}^{N} | Ω_{k} | \sqrt{∥ (\nabla δ x_{ROI}) |_{Ω_{k}} ∥^{2} + β},

is a differentiable approximation of the isotropic TV functional⁵⁵ and

{(\nabla δ x_{ROI}) |}_{Ω_{k}}

is the gradient of

δ x_{ROI}

at element

Ω_{k}

.

β > 0

is a small parameter that ensures TV

(δ x_{ROI})

is differentiable and

α

is the regularization parameter. In this paper,

α

and

β

were manually selected. For systematic approaches of the regularization parameter selection see, e.g., Refs. 56 and 57. The values of

α

and

β

used in our reconstructions are listed in Table 2.

Table 2

Regularization parameters of the TV regularization for δx. The first row shows α and β values of μa and μs′ used in the 2-D reconstructions. The second row shows the regularization parameters for the 3-D reconstructions.

	αδμa	αδμs	βδμa	βδμs
2-D	30,000	30	$5 \times 10^{- 5}$	$5 \times 10^{- 4}$
3-D	500	0.5	$2 \times 10^{- 6}$	$1 \times 10^{- 4}$

For modeling of $δ x$ in the conventional linear difference imaging Eq. (25), we used the same quadratic Gaussian smoothness regularization as before,

Eq. (28)

f (δ x) = ∥ L_{δ x} (δ x) ∥^{2},

except in this case, while modeling the random field

δ x

, the spatially constant part

x_{1, bg}

had a zero mean

x_{1, bg} = q I

with

q

\sim N (0, σ_{bg, x_{1}}^{2})

. The standard deviations of the background and inhomogeneities are given in Table 1.

Note that the standard deviations chosen for $x_{1}$ are larger than in $δ x$ since $x_{1}$ are absolute values of optical parameters unlike $δ x$ . Also, since $x_{1}$ is the “unchanging” parameter between observations $y_{1}$ and $y_{2}$ in the model (18), the presence of the same modeling errors in both observations $y_{1}$ and $y_{2}$ should mainly affect the estimate $x_{1}$ (not $δ x$ ). Thus, to allow for large variations in $x_{1}$ in the presence of modeling errors, we choose larger std’s for prior modeling of $x_{1}$ .

3.3.

Two-Dimensional Target and Simulations

The measurements $y_{1}$ and $y_{2}$ were simulated using 2-D simulation target states $x_{1}$ and $x_{2}$ shown in Fig. 1. The shape of the domain was extracted from a segmented adult brain CT scan. The diameter along the saggital plane was scaled to 100 mm (approximately the size of a newborn baby head). The state $x_{1}$ had three outer layers (mimicking the skin, skull, and cerebrospinal fluid) and two overlapping circular inclusions in the brain area. State $x_{2}$ was the same as $x_{1}$ , except that one additional inclusion in $μ_{a}$ and one in $μ_{s}^{'}$ were added in the dorsal (back) part of the brain. The optical properties of the states $x_{1}$ and $x_{2}$ are listed in Table 3. The measurement setup consisted of 16 sources and 16 detectors modeled as 1-mm-wide surface patches located on the boundary $\partial Ω$ . Random measurement noise $e_{i}, i = 1,2$ , drawn from zero-mean Gaussian distributions $π (e_{i}) = N (0, Γ_{e_{i}})$ , $i = 1,2$ where the standard deviations were specified as 1% of the simulated noise free measurement data, were added to the simulated measurement data. The means $e_{i, *} = 0$ , $i = 1,2$ and covariances $Γ_{e_{i}}$ , $i = 1,2$ were assumed known.

Fig. 1

Two-dimensional target optical properties: (a) absorption and (b) scattering. First column: reference state $x_{1}$ , second column: state $x_{2}$ , third column: change $δ x$ .

Table 3

Optical properties of the target head.

	μa(mm−1)	μs′(mm−1)
Layer 1	0.01	1
Layer 2	0.015	1.19
Layer 3	0.004	0.6
Inclusion 1	0.004	0.4
Inclusion 2	0.015	1.5
Inclusion 3 (only in $x_{2}$ )	0.018	1
Inclusion 4 (only in $x_{2}$ )	0.01	1.8
Background	0.01	1

3.3.1.

Using different optode arrangements on the boundary

We investigated situations where the optodes (16 sources, 16 detectors) were arranged at (a) equiangular intervals around the boundary of the 2-D target and at (b) equiangular intervals only in the dorsal part of the boundary of the 2-D target. In this case, we used a quadratic smoothness promoting functional for modeling $x_{1}$ , $f_{1} (x_{1}) = {∥ L_{x_{1}} (x_{1} - x_{1, *}) ∥}^{2}$ and sparsity promoting TV functional for modeling $δ x_{ROI}$ , $f_{2} (δ x_{ROI}) = α TV (δ x_{ROI})$ in the estimate (24). We used a quadratic smoothness promoting functional for modeling $δ x$ , $f (δ x) = {∥ L_{δ x} δ x ∥}^{2}$ in the estimate (25).

Figures 2(a) and 2(b) show the estimated optical coefficients from a simulation setup where the optodes were placed at equiangular intervals around the whole boundary. Panel (a) shows estimates of $x_{1}$ and $δ x_{ROI}$ using nonlinear difference imaging, Eq. (24). The ROI in the computations was selected as the dorsal hemisphere of the brain domain. The computation time ( $t_{CPU}$ ) of the nonlinear estimate was $t_{CPU} = 24.39 s$ . The reconstructions converged after three GN iterations. Panel (b) shows the corresponding conventional linear difference imaging estimate of $δ x$ , Eq. (25). The computation time of the linear estimate was $t_{CPU} = 1.31 s$ . The $t_{CPU}$ :s of all the remaining 2-D nonlinear and linear estimates were of the same magnitude as the values reported here. Figure 2(c) shows estimates $x_{1}$ and $δ x_{ROI}$ , Eq. (24) where the optodes were placed at the dorsal (upper) half of the boundary. Figure 2(d) shows the corresponding estimate $δ x$ using a linear difference imaging, Eq. (25). The error in the reconstructions

Error (\hat{x}) = \frac{{∥ \hat{x} - x_{0} ∥}^{2}}{{∥ x_{0} ∥}^{2}} \times 100 %,

where

\hat{x}

is estimated

δ μ_{a}

or

δ μ_{s}^{'}

and

x_{0} = x_{2} - x_{1}

is the simulated (true) target change, are listed in Table 4.

Fig. 2

Difference imaging using different optode arrangements on the boundary. (a) and (b) Estimated optical coefficients with nonlinear and linear difference imaging, respectively, with the optodes placed at equiangular intervals around the boundary. (c) and (d) Estimated optical coefficients with optodes placed at equiangular intervals only at the dorsal part of the boundary. In each row, the first column is $μ_{a, 1}$ , the second column is $δ μ_{a}$ , the third column is $μ_{s, 1}$ , and the fourth column is $δ μ_{s}^{'}$ . In the $δ x$ images, the black lines indicate the domain boundaries and the green circles indicate the position of the inclusions. The blank (white) areas in nonlinear estimates of ( $δ μ_{a}$ , $δ μ_{s}^{'}$ ) are regions outside region of interest (ROI).

Table 4

Reconstruction errors.

Section		Method	Error(δμa)	Error(δμs′)
3.3.1	Optodes around	Nonlinear	81	69
		Linear	87	58
	Optodes dorsal	Nonlinear	62	57
		Linear	92	71
3.3.2	No ROI	Nonlinear	79	70
	Quadratic regularization for $δ x_{ROI}$	Nonlinear	58	52
	No ROI	Linear	92	71
	ROI	Linear	85	63
3.3.3	Domain truncation	Nonlinear	73	67
		Linear	98	57
	Incorrect coupling	Nonlinear	70	50
		Linear	90	71

As can be seen from Fig. 2 and Table 4, nonlinear difference imaging shows better localization and recovery of the change compared to the conventional linear reconstruction for both optode arrangements. The nonlinear estimates of the change are not as spatially spread as the linear estimates. The effect is especially evident when only the upper part of the boundary is covered by sources and detectors. Since only the partial sensor coverage would be usually available in practical head imaging, the remaining 2-D simulations using the head domain were carried out using sources and detectors only at the dorsal part of the boundary.

3.3.2.

Using different region of interest constraints and regularizations

The purpose of this test case is to demonstrate that the improvement of the proposed approach over the conventional linearized reconstruction is not only because of (1) ROI constraint and (2) TV regularization for $δ x$ . The results are shown in Fig. 3. Figure 3(a) shows nonlinear difference reconstruction Eq. (24) without any ROI constraint i.e., $Ω_{ROI} = Ω$ . Figure 3(b) shows nonlinear difference reconstruction Eq. (24) using quadratic smoothness promoting functional for modeling $δ x_{ROI}$ , $f (δ x_{ROI}) = {∥ L_{δ x_{ROI}} δ x_{ROI} ∥}^{2}$ . Figure 3(c) shows conventional linear difference reconstruction Eq. (25), and Fig. 3(d) shows conventional linear difference reconstruction Eq. (25) with ROI constraint

\hat{δ x_{ROI}} = \arg \min_{δ x_{ROI}} {{∥ L_{δ e} (δ y - J M δ x_{ROI}) ∥}^{2} + f_{δ x_{ROI}} (δ x_{ROI})} .

Fig. 3

(a) Estimated optical coefficients using nonlinear difference imaging with no ROI constraint. (b) Estimated optical coefficients with nonlinear difference imaging using quadratic smoothness regularization for modeling $δ x_{ROI}$ . (c) Estimated change in optical coefficients using standard difference imaging. (d) Estimates obtained using ROI constant in standard difference imaging.

We can see that the results with nonlinear difference imaging [Figs. 3(a) and 3(b)] are quite similar to those obtained when using ROI constraint and TV regularization [Fig. 2(c)]. Also, the difference imaging estimates do not improve much by adding the ROI constraint. This result shows that the improvement of reconstruction is not only due to the choice of regularization or using the ROI constraint alone—it is in significant part due to the specific parametrization and formulation of the nonlinear difference imaging problem.

3.3.3.

Tolerance toward modeling errors

We tested with simulations the tolerance of the nonlinear and linear difference imaging toward modeling errors. Reconstructions in the presence of (a) domain truncation, (b) unknown optode coupling, and (c) unknown object shape were considered. The reconstructions are shown in Fig. 4. The errors in the reconstructions are listed in Table 4.

Fig. 4

Tolerance toward modeling errors. (a) and (b) Estimates of optical coefficients using the nonlinear difference imaging and the linear difference imaging in the presence of domain truncation. (c) and (d) Estimates in the presence of coupling errors. (e) and (f) Estimates in the presence of domain shape error.

Domain truncation: the first and second rows in Fig. 4 show estimates $x_{1}$ , $δ x_{ROI}$ and $δ x$ , Eqs. (24) and (25) in the presence of domain truncation, i.e., using a truncated domain as the model domain. The errors in the reconstructions are listed in Table 4.

Unknown optode coupling: the second and third rows in Fig. 4 show estimates $x_{1}$ , $δ x_{ROI}$ , and $δ x$ , Eq. (24) and (25) in the case of an inaccurately known optode coupling. The errors in the reconstructions are listed in Table 4. The coupling error was simulated as follows. Let $ζ = {(s, δ, d, η)}^{T} \in R^{512}$ represent the coupling coefficients of the 16 sources’ and 16 detectors’ amplitude and phases. Let us define a vector valued mapping $g (ζ) \in C^{256}$ such that

g_{k} (ζ) := {\hat{s}}_{p} {\hat{d}}_{q} = d_{p} s_{q} \exp [j (η_{p} + δ_{q})], k = (q - 1) 16 + p,

where

p

and

q

are the source and detector indexes, respectively. The coupling error

ϵ (ζ)

is given by²⁷

ϵ (ζ) = (\begin{matrix} Re \log [g (ζ)] \\ Im \log [g (ζ)] \end{matrix}) .

The data were corrupted with coupling error as $y_{i, ϵ (ζ)} = y_{i} + ϵ (ζ)$ , $i = 1,2$ . In this case, the source and detector amplitude and phase coefficients were drawn from uniform distribution $s_{p}$ , $d_{q} \sim U (0.9, 1)$ , $δ_{p}$ , $η_{q} \sim U (0, π / 360)$ . The data used in estimate (24) was $\tilde{y} = {(y_{1, ϵ (ζ)}, y_{2, ϵ (ζ)})}^{T}$ and for estimate (25), it was $δ y = y_{2, ϵ (ζ)} - y_{1, ϵ (ζ)}$ .

Unknown object shape: the fourth and fifth rows in Fig. 4 show estimates $x_{1}$ , $δ x_{ROI}$ , and $δ x$ , Eqs. (24) and (25) in the case of using incorrect model domain, i.e., in this case an incorrectly shaped domain (domain obtained using some other segmented adult CT scan) was used as the model domain. The reconstruction errors for this case are not listed in Table 4 since the deformation map of the true optical properties from the measurement domain to the model domain is not known.

From the estimates shown in Fig. 4, we can observe that in most cases, the linear difference reconstruction Eq. (25) is indicative of the location of change, although the spatial resolution is weak. In the estimates obtained using the nonlinear difference imaging Eq. (24), the reconstructed initial state $x_{1}$ is heavily affected by the modeling errors; however, the estimates $δ x$ are relatively free from artifacts. The reason for the modeling errors not affecting significantly the estimates $δ x$ lies in the parametrization of $x_{2}$ as a linear combination of the initial state $x_{1}$ and the change $δ x$ as the unknowns. The variable $x_{1}$ is invariant between models $y_{1} = A (x_{1}) + e_{1}$ and $y_{2} = A (x_{1} + M δ x_{ROI}) + e_{2}$ . Consequently, when the proposed parametrization is used, the errors caused by the (invariant) modeling errors are propagated mainly to the estimate of $x_{1}$ , which consists the parameters that are common for the models of both observations $y_{1}$ and $y_{2}$ .

3.4.

Reconstructions in Three Dimensions with Experimental Data

3.4.1.

Experimental details

The experiment was carried out with the frequency domain DOT instrument at the Aalto University, Helsinki.⁵⁸ The measurement domains were cylinders with radius 35 mm and height 110 mm (see Fig. 5). The cylindrical phantoms corresponding to states $x_{1}$ and $x_{2}$ are illustrated in Fig. 5. The background optical properties were approximately $μ_{a, bg} = 0.01 {mm}^{- 1}$ and $μ_{s, bg}^{'} = 1 {mm}^{- 1}$ at wavelength 800 nm for both phantoms. The cylindrical inclusions in $x_{2}$ , which both have diameter and height of 9.5 mm, are located such that the central planes of the inclusions coincide with the central $x y$ -plane of the cylinder domain. The optical properties of inclusion 1 are approximately $μ_{a, inc . 1} = 0.02 {mm}^{- 1}$ , $μ_{s, inc . 1}^{'} = 1 {mm}^{- 1}$ (i.e., purely absorption contrast) and the properties of inclusion 2 are $μ_{a, inc . 2} = 0.01 {mm}^{- 1}$ , $μ_{s, inc . 2}^{'} = 2 {mm}^{- 1}$ (i.e., purely scatter contrast), respectively. Absolute imaging reconstructions using the phantom with inclusions are presented in Refs. 53, 58, and 59.

Fig. 5

(a) Photograph of the diffuse optical tomography experiment using one of the phantoms. (b) Cylindrical phantoms $x_{1}$ and $x_{2}$ with the position of sources (red circles) and detectors (blue crosses). $x_{1}$ is a homogeneous reference phantom with optical coefficients approximately $μ_{a, bg} = 0.01 {mm}^{- 1}$ and $μ_{s, bg}^{'} = 1 {mm}^{- 1}$ at wavelength 800 nm. $x_{2}$ has the same background optical coefficients as $x_{1}$ and it has two additional inclusion with optical coefficients approximately $μ_{a, inc . 1} = 0.02 {mm}^{- 1}$ , $μ_{s, inc . 1}^{'} = 1 {mm}^{- 1}$ (i.e., purely absorption contrast) of inclusion 1 and $μ_{a, i nc . 2} = 0.01 {mm}^{- 1}$ , $μ_{s, inc . 2}^{'} = 2 {mm}^{- 1}$ (i.e., purely scatter contrast) of inclusion 2. The gray patches on the cylinders show the ROI.

The source and detector configuration in the experiment consisted of 16 sources and 15 detectors arranged in interleaved order on two rings located 6 mm above and below the central $x y$ -plane of the cylinder domain. The locations of sources and detectors are shown with red and blue circles respectively in Fig. 5. The measurements were carried out at 785 nm with an optical power of 8 mW and modulation frequency $ω = 100 MHz / (2 π)$ . The log amplitude and phase shift of the transmitted light was recorded and the nearest measurement data from each source position were removed, leading to real valued measurement vectors $y_{i} \in R^{360}$ , $i = 1,2$ . We employed an error model

e \sim N (0, Γ_{e_{i}}), i = 1,2,

where the square roots of the diagonal elements (standard deviations) of error covariances

Γ_{e_{i}}

were specified as 1% of the absolute values of

y_{i}

,

i = 1,2

, implying that the standard deviations of the measurement errors are assumed to be 1% of the measured absolute values of the log(amplitude) and phase.

3.4.2.

Data calibration and initialization of the estimation

Raw instrument data cannot be directly used as a direct equivalent to simulated data from a model. The measured phase and amplitude were calibrated using the procedure described in Ref. 58, accounting for differences between the lengths and coupling efficiencies between different source and detector channels, as well as for the effects of detector gain adjustment during the data collection. Finally, to calibrate the forward model to the measurement setup, the initial estimates for the optical properties of the model were assumed to be homogeneous and the measured data was used to fit global values of absorption and scattering in the model as well as global coupling factors for phase and amplitude between model and measurement. The coupling is intrinsically unknown a priori; therefore, the measured data are corrected by the coupling coefficient obtained from this optimization process. In a sense, the instrument data are calibrated to match the absolute values of the simulations from the forward model. However, to avoid unrealistical calibration setup of having access to data from the homogeneous target, the four-parameter calibration was carried out using the data $y_{2}$ , corresponding to the nonhomogeneous state after the change. Following the initial estimation procedure in Ref. 59, the log of source strength and phase coupling between the modeled phase and log amplitude, and the measured phase and log amplitude were modeled by additive constants. In other words, we assumed that the coupling factors are constant for all source and detector fibers. Thus, the initialization step consisted of a four-parameter fit of global background parameters $μ_{a, 0}$ and $μ_{s, 0}^{'}$ as well as a global additive shift $η$ of log amplitude data and a global additive shift $φ$ of phase data

{μ_{a, 0}, μ_{s, 0}^{'}, η, φ} = \arg \min_{μ_{a, 0}, μ_{s, 0}^{'}, η, φ} {{∥ L_{e_{2}} [(y_{2} + Δ y) - A (μ_{a, 0}, μ_{s, 0}^{'})] ∥}^{2}},

where

Δ y = {(η, φ)}^{T}

. The initialization problem was solved by a GN optimization method with an explicit line search.⁵³ Once the initialization was completed, the measurement data were transformed for the nonlinear estimation Eq. (24) by the recovered global offsets as

\tilde{y} = {(y_{1} + Δ y, y_{2} + Δ y)}^{T}

, and the initial parameter values for the nonlinear estimation were set to the estimated values

μ_{a, 0}

and

μ_{s, 0}^{'}

.

3.4.3.

Results

Figure 6 shows 2-D slices of 3-D reconstructions obtained using nonlinear difference imaging Eq. (24). The ROI in this case was selected as a central part of height $z = 22 mm$ of the cylinder (see Fig. 5). Figs. 6(a) and 6(b) are the absorption images and Figs. 6(c) and 6(d) are the scattering images. The computation time of the nonlinear estimate was $t_{CPU} = 3247.68 s$ . Figure 7 shows the corresponding 3-D reconstructions with linear difference imaging Eq. (25). The computation time of the linear estimate was $t_{CPU} = 87.95 s$ .

Fig. 6

Estimates using nonlinear difference imaging: (a), (b) $μ_{a}$ and (c), (d) $μ_{s}^{'}$ . (a) and (c) The horizontal slices of the phantom along $z$ -axis at $z = 22, 55, 99 mm$ . (b) Vertical slices along $x$ -axis at $x = 14, 35, 63 mm$ . (d) Vertical slices along $y$ -axis at $y = 14, 35, 63 mm$ . The blank (white) areas are regions outside ROI.

Fig. 7

Estimates using linear difference imaging. (a), (b) $μ_{a}$ and (c), (d) $μ_{s}^{'}$ . (a) and (c) Horizontal slices of the phantom (height 110 mm, diameter 70 mm) along $z$ -axis at $z = 22, 55, 99 mm$ . (b) Vertical slices along $x$ -axis at $x = 14, 35, 63 mm$ . (d) Vertical slices along $y$ -axis at $y = 14, 35, 63 mm$ .

From Figs. 6 and 7, one can see that the proposed approach of nonlinear difference imaging shows better recovery of the inclusions when compared to the conventional linear difference imaging. Also, there is no “ringing” effect in nonlinear reconstruction, while there are severe ones in linear reconstruction. These could result in a masking effect when more target inclusions are present. These results are in agreement with the 2-D simulation results.

In the 3-D case, the CPU time needed for the reconstruction with the nonlinear approach was approximately 36 times the time needed for the linear reconstruction. This can be potentially problematic when processing long time series of data, especially if the reconstructions are needed online. However, in many cases the requirement of better reconstruction accuracy may outweigh the disadvantage of longer CPU time, especially if it is sufficient to obtain the results off-line. It should be also noted that for both the linear and nonlinear approach the reported CPU times are based on nonoptimized implementation on MATLAB® using the TOAST toolbox. Where needed, the computations can be made faster by optimizing the implementation.

4. Conclusions

We applied a new approach to difference imaging in DOT. In the approach, the optical coefficients after a change are parameterized as a linear combination of the (unknown) initial state and the change in optical properties. The DOT measurements taken before and after the change are concatenated into a single measurement vector and the inverse problem is stated as finding the initial optical coefficients and the change given the combined data. This model allows the use of separate spatial models for the initial state and the change in optical coefficients. Furthermore, it allows the use of a ROI constraint for the change in optical coefficients in a straightforward way. The approach was also tested with 2-D simulations using different optode placement settings and also in the presence of modeling errors arising from domain truncation, unknown optode coupling and unknown domain shape. The conventional linearized difference imaging was used as a reference approach. The approach was tested with experimental phantom data. The results show that the proposed approach produces better reconstructions compared to standard linear difference imaging and the approach is robust for the modeling errors at least in a similar extent as the conventional linear reconstruction approach. We believe that the proposed approach can improve the accuracy of difference imaging compared to the linearization-based approach, which is the current standard in difference imaging.

Acknowledgments

This work was supported by the Academy of Finland (projects 136220, 272803, 269282, and 250215 Finnish Centre of Excellence in Inverse Problems Research) and Finnish Doctoral Programme in Computational Sciences. The authors would also like to acknowledge Dong Liu for useful discussions.

References

1.

J. P. Culver et al., “Volumetric diffuse optical tomography of brain activity,” Opt. Lett., 28 (21), 2061 –2063 (2003). http://dx.doi.org/10.1364/OL.28.002061 OPLEDP 0146-9592 Google Scholar

2.

J. C. Hebden et al., “Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol., 47 (23), 4155 –4166 (2002). http://dx.doi.org/10.1088/0031-9155/47/23/303 PHMBA7 0031-9155 Google Scholar

3.

T. Näsi et al., “Effect of task-related extracerebral circulation on diffuse optical tomography: experimental data and simulations on the forehead,” Biomed. Opt. Express, 4 (3), 412 –426 (2013). http://dx.doi.org/10.1364/BOE.4.000412 BOEICL 2156-7085 Google Scholar

4.

B. J. Tromberg et al., “Assessing the future of diffuse optical imaging technologies for breast cancer management,” Med. Phys., 35 (6), 2443 (2008). http://dx.doi.org/10.1118/1.2919078 MPHYA6 0094-2405 Google Scholar

5.

D. Leff et 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

6.

G. Gulsen et al., “Combined diffuse optical tomography (DOT) and MRI system for cancer imaging in small animals,” Technol. Cancer Res. Treat., 5 (4), 351 –363 (2006). http://dx.doi.org/10.1177/153303460600500407 Google Scholar

7.

J. P. Culver et al., “Diffuse optical tomography of cerebral blood flow, oxygenation and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab., 23 (8), 911 –924 (2003). http://dx.doi.org/10.1097/01.WCB.0000076703.71231.BB Google Scholar

8.

S. R. Hintz et al., “Bedside functional imaging of the premature infant brain during passive motor activation,” J. Perinat. Med., 29 (4), 335 –343 (2001). http://dx.doi.org/10.1515/JPM.2001.048 Google Scholar

9.

A. Gibson et al., “Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate,” NeuroImage, 30 (2), 521 –528 (2006). http://dx.doi.org/10.1016/j.neuroimage.2005.08.059 NEIMEF 1053-8119 Google Scholar

10.

J. C. Hebden et al., “Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys. Med. Biol., 49 (7), 1117 (2004). http://dx.doi.org/10.1088/0031-9155/49/7/003 PHMBA7 0031-9155 Google Scholar

11.

T. Zhang et al., “Pre-seizure state identified by diffuse optical tomography,” Sci. Rep., 4 3798 (2014). Google Scholar

12.

A. T. Eggebrecht et al., “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” NeuroImage, 61 (4), 1120 –1128 (2012). http://dx.doi.org/10.1016/j.neuroimage.2012.01.124 NEIMEF 1053-8119 Google Scholar

13.

Y. Zhan et al., “Image quality analysis of high-density diffuse optical tomography incorporating a subject-specific head model,” Front. Neuroenerg., 4 (6), (2012). http://dx.doi.org/10.3389/fnene.2012.00006 Google Scholar

14.

E. M. C. Hillman et al., “Calibration techniques and datatype extraction for time-resolved optical tomography,” Rev. Sci. Instrum., 71 (9), 3415 (2000). http://dx.doi.org/10.1063/1.1287748 RSINAK 0034-6748 Google Scholar

15.

C. H. Schmitz et al., “Instrumentation and calibration protocol for imaging dynamic features in dense-scattering media by optical tomography,” Appl. Opt., 39 (34), 6466 –6486 (2000). http://dx.doi.org/10.1364/AO.39.006466 APOPAI 0003-6935 Google Scholar

16.

D. Boas, T. Gaudette and S. Arridge, “Simultaneous imaging and optode calibration with diffuse optical tomography,” Opt. Express, 8 (5), 263 –270 (2001). http://dx.doi.org/10.1364/OE.8.000263 OPEXFF 1094-4087 Google Scholar

17.

M. Schweiger et al., “Image reconstruction in optical tomography in the presence of coupling errors,” Appl. Opt., 46 (14), 2743 –2756 (2007). http://dx.doi.org/10.1364/AO.46.002743 APOPAI 0003-6935 Google Scholar

18.

X. Wu et al., “Quantitative evaluation of atlas-based high-density diffuse optical tomography for imaging of the human visual cortex,” Biomed. Opt. Express, 5 (11), 3882 –3900 (2014). http://dx.doi.org/10.1364/BOE.5.003882 BOEICL 2156-7085 Google Scholar

19.

A. Gibson et al., “Optical tomography of a realistic neonatal head phantom,” Appl. Opt., 42 (16), 3109 –3116 (2003). http://dx.doi.org/10.1364/AO.42.003109 APOPAI 0003-6935 Google Scholar

20.

A. H. Barnett et al., “Robust inference of baseline optical properties of the human head with three-dimensional segmentation from magnetic resonance imaging,” Appl. Opt., 42 (16), 3095 –3108 (2003). http://dx.doi.org/10.1364/AO.42.003095 APOPAI 0003-6935 Google Scholar

21.

B. W. Pogue and K. D. Paulsen, “High-resolution near-infrared tomographic imaging simulations of the rat cranium by use of a priori magnetic resonance imaging structural information,” Opt. Lett., 23 (21), 1716 –1718 (1998). http://dx.doi.org/10.1364/OL.23.001716 OPLEDP 0146-9592 Google Scholar

22.

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

23.

S. R. Arridge et al., “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl., 22 (1), 175 –195 (2006). http://dx.doi.org/10.1088/0266-5611/22/1/010 Google Scholar

24.

T. Tarvainen et al., “An approximation error approach for compensating for modelling errors between the radiative transfer equation and the diffusion approximation in diffuse optical tomography,” Inverse Probl., 26 (1), 015005 (2010). http://dx.doi.org/10.1088/0266-5611/26/1/015005 Google Scholar

25.

T. Tarvainen et al., “Corrections to linear methods for diffuse optical tomography using approximation error modelling,” Biomed. Opt. Express, 1 (1), 209 –222 (2010). http://dx.doi.org/10.1364/BOE.1.000209 BOEICL 2156-7085 Google Scholar

26.

V. Kolehmainen et al., “Marginalization of uninteresting distributed parameters in inverse problems–application to diffuse optical tomography,” Int. J. Uncertainty Quantif., 1 (1), 1 –17 (2011). http://dx.doi.org/10.1615/Int.J.UncertaintyQuantification.v1.i1.10 Google Scholar

27.

M. Mozumder et al., “Compensation of optode sensitivity and position errors in diffuse optical tomography using the approximation error approach,” Biomed. Opt. Express, 4 (10), 2015 –2031 (2013). http://dx.doi.org/10.1364/BOE.4.002015 BOEICL 2156-7085 Google Scholar

28.

M. Mozumder et al., “Compensation of modeling errors due to unknown domain boundary in diffuse optical tomography,” J. Opt. Soc. Am. A, 31 (8), 1847 –1855 (2014). http://dx.doi.org/10.1364/JOSAA.31.001847 JOAOD6 0740-3232 Google Scholar

29.

Y. Pei, F. B. Lin and R. Barbour, “Modeling of sensitivity and resolution to an included object in homogeneous scattering media and in MRI-derived breast maps,” Opt. Express, 5 (10), 203 –219 (1999). http://dx.doi.org/10.1364/OE.5.000203 OPEXFF 1094-4087 Google Scholar

30.

A. Bluestone et 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

31.

Y. Pei, H. L. Graber and R. L. Barbour, “Influence of systematic errors in reference states on image quality and on stability of derived information for DC optical imaging,” Appl. Opt., 40 (31), 5755 –5769 (2001). http://dx.doi.org/10.1364/AO.40.005755 APOPAI 0003-6935 Google Scholar

32.

T. Austin et al., “Three dimensional optical imaging of blood volume and oxygenation in the neonatal brain,” NeuroImage, 31 (4), 1426 –1433 (2006). http://dx.doi.org/10.1016/j.neuroimage.2006.02.038 NEIMEF 1053-8119 Google Scholar

33.

R. L. Barbour et al., “Optical tomographic imaging of dynamic features of dense-scattering media,” J. Opt. Soc. Am. A, 18 (12), 3018 –3036 (2001). http://dx.doi.org/10.1364/JOSAA.18.003018 JOAOD6 0740-3232 Google Scholar

34.

Y. Xu, H. L. Graber and R. L. Barbour, “Image correction algorithm for functional three-dimensional diffuse optical tomography brain imaging,” Appl. Opt., 46 (10), 1693 –1704 (2007). http://dx.doi.org/10.1364/AO.46.001693 APOPAI 0003-6935 Google Scholar

35.

R. L. Barbour et al., “Strategies for imaging diffusing media,” Transp. Theory Stat. Phys., 33 (3–4), 361 –371 (2004). http://dx.doi.org/10.1081/TT-200051950 TTSPB4 0041-1450 Google Scholar

36.

H. L. Graber et al., “Spatial deconvolution technique to improve the accuracy of reconstructed three-dimensional diffuse optical tomographic images,” Appl. Opt., 44 (6), 941 –953 (2005). http://dx.doi.org/10.1364/AO.44.000941 APOPAI 0003-6935 Google Scholar

37.

Y. Xu et al., “Improved accuracy of reconstructed diffuse optical tomographic images by means of spatial deconvolution: two-dimensional quantitative characterization,” Appl. Opt., 44 (11), 2115 –2139 (2005). http://dx.doi.org/10.1364/AO.44.002115 APOPAI 0003-6935 Google Scholar

38.

Y. Xu et al., “Image quality improvement via spatial deconvolution in optical tomography: time-series imaging,” J. Biomed. Opt., 10 (5), 051701 (2005). http://dx.doi.org/10.1117/1.2103747 JBOPFO 1083-3668 Google Scholar

39.

H. L. Graber, Y. Xu and R. L. Barbour, “Image correction scheme applied to functional diffuse optical tomography scattering images,” Appl. Opt., 46 (10), 1705 –1716 (2007). http://dx.doi.org/10.1364/AO.46.001705 APOPAI 0003-6935 Google Scholar

40.

J. Heiskala, P. Hiltunen and I. Nissilä, “Significance of background optical properties, time-resolved information and optode arrangement in diffuse optical imaging of term neonates,” Phys. Med. Biol., 54 (3), 535 (2009). http://dx.doi.org/10.1088/0031-9155/54/3/005 PHMBA7 0031-9155 Google Scholar

41.

D. Liu et al., “A nonlinear approach to difference imaging in EIT; assessment of the robustness in the presence of modelling errors,” Inverse Probl., 31 (3), 035012 (2015). http://dx.doi.org/10.1088/0266-5611/31/3/035012 INPEEY 0266-5611 Google Scholar

42.

D. Liu et al., “Estimation of conductivity changes in a region of interest with electrical impedance tomography,” Inverse Probl. Imaging, 9 (1), 211 –229 (2015). http://dx.doi.org/10.3934/ipi.2015.9.211 Google Scholar

43.

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

44.

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

45.

A. Douiri et al., “Anisotropic diffusion regularization methods for diffuse optical tomography using edge prior information,” Meas. Sci. Technol., 18 (1), 87 (2007). http://dx.doi.org/10.1088/0957-0233/18/1/011 MSTCEP 0957-0233 Google Scholar

46.

P. K. Yalavarthy et al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express, 15 (13), 8043 –8058 (2007). http://dx.doi.org/10.1364/OE.15.008043 OPEXFF 1094-4087 Google Scholar

47.

K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization,” Appl. Opt., 35 (19), 3447 –3458 (1996). http://dx.doi.org/10.1364/AO.35.003447 APOPAI 0003-6935 Google Scholar

48.

V. Kolehmainen, “Novel approaches to image reconstruction in diffusion tomography,” University of Kuopio, Kuopio, Finland, (2001). Google Scholar

49.

N. Cao, A. Nehorai and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express, 15 (21), 13695 –13708 (2007). http://dx.doi.org/10.1364/OE.15.013695 OPEXFF 1094-4087 Google Scholar

50.

J. F. P.-J. Abascal et al., “Fluorescence diffuse optical tomography using the split Bregman method,” Med. Phys., 38 (11), 6275 –6284 (2011). http://dx.doi.org/10.1118/1.3656063 MPHYA6 0094-2405 Google Scholar

51.

J. C. Ye et al., “Optical diffusion tomography by iterative-coordinate-descent optimization in a Bayesian framework,” J. Opt. Soc. Am. A, 16 (10), 2400 –2412 (1999). http://dx.doi.org/10.1364/JOSAA.16.002400 JOAOD6 0740-3232 Google Scholar

52.

D. A. Boas, A. M. Dale and M. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” NeuroImage, 23 (Suppl 1), S275 –S288 (2004). http://dx.doi.org/10.1016/j.neuroimage.2004.07.011 NEIMEF 1053-8119 Google Scholar

53.

M. Schweiger, S. R. Arridge and I. Nissilä, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol., 50 (10), 2365 –2386 (2005). http://dx.doi.org/10.1088/0031-9155/50/10/013 PHMBA7 0031-9155 Google Scholar

54.

J. Heiskala et al., “Approximation error method can reduce artifacts due to scalp blood flow in optical brain activation imaging,” J. Biomed. Opt., 17 (9), 096012 (2012). http://dx.doi.org/10.1117/1.JBO.17.9.096012 JBOPFO 1083-3668 Google Scholar

55.

L. I. Rudin, S. Osher and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D Nonlinear Phenom., 60 (1–4), 259 –268 (1992). http://dx.doi.org/10.1016/0167-2789(92)90242-F Google Scholar

56.

Y. Mamatjan et al., “An experimental clinical evaluation of EIT imaging with l-1 data and image norms,” Physiol. Meas., 34 (9), 1027 (2013). http://dx.doi.org/10.1088/0967-3334/34/9/1027 PMEAE3 0967-3334 Google Scholar

57.

A. S. K. Karhunen et al., “Electrical resistance tomography for assessment of cracks in concrete,” ACI Mater. J., 107 (5), 523 –531 (2010). Google Scholar

58.

I. Nissilä et al., “Instrumentation and calibration methods for the multichannel measurement of phase and amplitude in optical tomography,” Rev. Sci. Instrum., 76 (4), 044302 (2005). http://dx.doi.org/10.1063/1.1884193 RSINAK 0034-6748 Google Scholar

59.

V. Kolehmainen et al., “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A., 26 (10), 2257 –2268 (2009). http://dx.doi.org/10.1364/JOSAA.26.002257 JOAOD6 0740-3232 Google Scholar

Biography

Meghdoot Mozumder is a doctoral student at the University of Eastern Finland. He received his MSc degree in physics from the Indian Institute of Technology, Kanpur, in 2011. His current research interests include optical imaging and inverse problems.

Tanja Tarvainen is an academy research fellow at the University of Eastern Finland. She received her PhD in 2006 from the University of Kuopio. She is the author of more than 30 journal papers. Her current research interests include diffuse optical tomography, quantitative photoacoustic tomography, and Bayesian inverse problems.

Aku Seppänen is an academy research fellow at the University of Eastern Finland. He received his PhD degree in 2006 from the University of Kuopio. He has authored 27 journal papers and two book chapters. His research focuses on inverse problems, with applications to imaging, nondestructive testing, and remote sensing.

Ilkka Nissilä is an Academy Research Fellow at Aalto University. He received his DSc degree at Helsinki University of Technology in 2004. His main research interests includes the development of diffuse optical tomography in combination with other modalities for the neuroimaging of children.

Simon R. Arridge is a professor of image processing at the Department of Computer Science and the director for the Centre for Inverse Problems at University College London. He received his PhD in 1990 from University College London. He is an author on more than 150 journal papers. His research interests include inverse problems methods, image reconstruction, regularization, modeling, and numerical methods.

Ville Kolehmainen is an associate professor at the University of Eastern Finland. He received his PhD in 2001 from the University of Kuopio. He is the author of more than 50 journal papers and has written four book chapters. His current research interests include computational inverse problems, especially in tomography imaging.

Citation Download Citation

Meghdoot Mozumder, Tanja Tarvainen, Aku Seppänen, Ilkka Nissilä, Simon R. Arridge, and Ville Kolehmainen "Nonlinear approach to difference imaging in diffuse optical tomography," Journal of Biomedical Optics 20(10), 105001 (6 October 2015). https://doi.org/10.1117/1.JBO.20.10.105001

Published: 6 October 2015

Access the abstract

JOURNAL ARTICLE
12 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 10 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Nonlinear optics

Data modeling

Error analysis

Diffuse optical imaging

Diffuse optical tomography

Sensors

Optical imaging

1.

Introduction

2.

Forward Model of Diffuse Optical Tomography

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

2.1.

Absolute Imaging

Eq. (9)

2.2.

Difference Imaging

Eq. (10)

Eq. (11)

2.2.1.

Conventional linear approach to difference imaging

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

2.2.2.

Nonlinear approach to difference imaging

Eq. (16)

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

Eq. (22)

Eq. (23)

3.

Results

3.1.

Estimates

Eq. (24)

Eq. (25)

3.2.

Regularization Functionals

Eq. (26)

Table 1

Eq. (27)

Table 2

Eq. (28)

3.3.

Two-Dimensional Target and Simulations

Fig. 1

Table 3

3.3.1.

Using different optode arrangements on the boundary

Fig. 2

Table 4

3.3.2.

Using different region of interest constraints and regularizations

Fig. 3

3.3.3.

Tolerance toward modeling errors

Fig. 4

3.4.

Reconstructions in Three Dimensions with Experimental Data

3.4.1.

Experimental details

Fig. 5

3.4.2.

Data calibration and initialization of the estimation

3.4.3.

Results

Fig. 6

Fig. 7

4.

Conclusions

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases