appliedsciences
ArticleA Hybrid Inversion Scheme Combining Markov
Chain Monte Carlo and Iterative Methods for
Determining Optical Properties of Random MediaYu Jiang 1 , Yoko Hoshi 2 , Manabu Machida 2,∗
1
2
3

*

and Gen Nakamura 3

School of Mathematics, Shanghai University of Finance and Economics, Shanghai 200433, China
Institute for Medical Photonics Research, Hamamatsu University School of Medicine,
Hamamatsu 431-3192, Japan
Department of Mathematics, Hokkaido University, Sapporo 060-0810, Japan
Correspondence: machida@hama-med.ac.jp


Received: 31 July 2019; Accepted: 21 August 2019; Published: 24 August 2019

Abstract: Near-infrared spectroscopy (NIRS) including diffuse optical tomography is an imaging
modality which makes use of diffuse light propagation in random media. When optical propertiesof a random medium are investigated from boundary measurements of reflected or transmitted
light, iterative inversion schemes such as the Levenberg–Marquardt algorithm are known to failwhen initial guesses are not close enough to the true value of the coefficient to be reconstructed.
In this paper, we investigate how this weakness of iterative schemes is overcome using Markov chain
Monte Carlo. Using time-resolved measurements performed against a polyurethane-based phantom,
we present a case that the Levenberg–Marquardt algorithm fails to work but the proposed hybrid
method works well. Then, with a toy model of diffuse optical tomography we illustrate that theLevenberg–Marquardt method fails when it is trapped by a local minimum but the hybrid method
can escape from local minima by using the Metropolis–Hastings Markov chain Monte Carlo algorithm
until it reaches the valley of the global minimum. The proposed hybrid scheme can be applied todifferent inverse problems in NIRS which are solved iteratively. We find that for both numerical andphantom experiments, optical properties such as the absorption and reduced scattering coefficientscan be retrieved without being trapped by a local minimum when Monte Carlo simulation is run
only about 100 steps before switching to an iterative method. The hybrid method is compared withsimulated annealing. Although the Metropolis–Hastings MCMC arrives at the steady state at about
10,000 Monte Carlo steps, in the hybrid method the Monte Carlo simulation can be stopped way
before the burn-in time.
Keywords: near-infrared spectroscopy; diffuse light; inverse problems; optical tomography

1. Introduction
In near-infrared spectroscopy (NIRS), we estimate optical properties of biological tissue by solving
inverse diffuse problems [1,2]. Such inverse problems are commonly solved by means of iterative
methods. In the case of a homogeneous medium, absorption and reduced scattering coefficients ofthe medium can be obtained with iterative methods such as the Levenberg–Marquardt algorithm[3,4] from time-resolved measurements (for example, see the review article [5]). Neuroimaging forthe human brain by NIRS through the neurovascular coupling has been developed and is called
functional NIRS (fNIRS) [6–9]. Since the region of interest on the head can be small, it is possible toassume a simple geometry of the half space. Quantitative measurements of inter-regional differencesin neuronal activity requires accurate estimates of optical properties.

Appl. Sci. 2019, 9, 3500; doi:10.3390/app9173500

www.mdpi.com/journal/applsci

Appl. Sci. 2019, 9, 3500

2 of 17

In Ref. [10], the Nelder–Mead simplex method was used to retrieve optical parameters in layered
tissue. Heterogeneity of optical properties can be obtained by diffuse optical tomography [1,11,12].
Iterative numerical schemes are used to minimize the cost function when solving these inverse
problems. A gradient-based approach [13] was used to detect breast cancer [14]. The brain activityof a newborn infant was investigated [15] with diffuse optical tomography by TOAST (temporal
optical absorption and scattering tomography) [16,17], in which iterative algorithms such as thenonlinear conjugate gradients, damped Gauss–Newton method, and Levenberg–Marquardt method
are implemented. Diffuse optical tomography was performed on human lower legs and a forearm withthe algebraic reconstruction algorithm in the framework of the modified generalized pulse spectrum
technique [18]. See Refs. [19,20] for numerical techniques for diffuse optical tomography. For these
iterative numerical schemes to work, it is important to choose a good initial guess for the initial valueof the iteration.
Solving inverse problems by the Bayesian approach has been sought as an alternative way.
In Ref. [21], a novel use of the Bayesian approach was considered to take modeling error into
account. The Bayesian inversion with the Metropolis–Hastings Markov chain Monte Carlo wasused in Refs. [22,23]. The Bayesian approach was used to determine optical parameters of the humanhead [24]. Although the Markov chain Monte Carlo (MCMC) approach is in principle able to escapefrom local minima, it is computationally time consuming. Hence, despite the above-mentioned efforts,
the use of Monte Carlo for inverse problems in NIRS has been extremely limited.
In this paper, we shed light on Markov chain Monte Carlo once again by combining it withan iterative scheme, and investigate the use of it in NIRS. In particular, we test a hybrid numerical
scheme of Markov chain Monte Carlo and an iterative method. In the proposed hybrid scheme,
the Markov chain Monte Carlo algorithm is first used to provide an initial guess using jumps in thelandscape of the cost function, and then an iterative method is used after the initial large fluctuation.
Thus, the proposed method realizes a fast reconstruction while, in the beginning, obtained values
at Monte Carlo steps jump around in the landscape of the cost function. The MCMC simulationis necessary only for the first 100 steps. Then the hybrid method starts to use an iterative scheme
and reconstructs optical properties by searching the global minimum. The computation time of theiterative scheme is negligible compared with that of the Monte Carlo simulation. Since the Monte
Carlo simulation can be stopped way before the burn-in time for the hybrid method, the proposed
scheme is at least ten times faster than simulated annealing, which is implemented by the naive use of
the Metropolis–Hastings Markov chain Monte Carlo.
The rest of the paper is organized as follows. We develop diffusion theory in Section 2.
A polyurethane-based phantom and numerical phantom are also described in Section 2. Section 3 isdevoted to experimental and numerical results. Finally, discussion is given in Section 4.
2. Materials and Methods
2.1. Diffusion Theory
2.1.1. Diffuse Light in Three Dimensions
Let us suppose that a random medium occupies the three-dimensional half space. We assume
that the random medium is characterized by the diffusion coefficient D and absorption coefficientµ a . We have D = 1/(3µ0s ), where µ0s is the reduced scattering coefficient. Position in the half space(−∞ < x < ∞, −∞ < y < ∞, 0 < z < ∞) is denoted by r = (ρ, z), ρ = ( x, y). Let t denote time. Let cbe the speed of light in the medium. We assume an incident beam on the boundary at the origin in thex-y plane. The energy density u of light in the medium obeys the following diffusion equation.

1 ∂− D∆ + µ a u(r, t) = 0,
c ∂t

(1)

Appl. Sci. 2019, 9, 3500

3 of 17

with the boundary condition−`∂u(r, t) + u(r, t) = δ( x )δ(y)q(t),∂z(2)

and the initial condition u(r, 0) = 0. The right-hand side of the boundary condition is the incident
laser beam which illuminates the medium at the origin (0, 0, 0) with the temporal profile q(t). In thephantom experiment described below, the incident light enters the phantom in the positive z direction.
The extrapolation distance ` is a nonnegative constant. If we consider the diffuse surface reflection, we
have [25]` = 2D

1 + rd,1 − rd(3)

where the internal reflection rd is given by [26] rd = −1.4399n−2 + 0.7099n−1 + 0.6681 + 0.0636n.
Let us consider the corresponding surface Green’s function Gs (r, t; ρ0 , s), which satisfies Equation (1)and the boundary condition


−`


∂
+ 1 Gs (r, t; ρ0 , s) = δ( x − x 0 )δ(y − y0 )δ(t − s),∂z(4)together with the initial condition Gs (r, 0; ρ0 , s) = 0. We obtain [27–30]"
2( x − x 0 )2 +(y−y0 )22e−µa c(t−s)
cD−− 4Dcz(t−s)
4Dc(t−s)eeGs (r, t; ρ , s) =` (4πDc(t − s))3/2
0

0 2

0 2

) +(y−y )z Dc(t−s)e−µa c(t−s) − (x−x4Dc
(t−s)
−
e
e ` e `2 erfc
4π ` Dc(t − s)z + 2Dc(t − s)/`p
4Dc(t − s)

(5)

!#
,

t > s,

with Gs (r, t; ρ0 , s) = 0 for t ≤ s. The complementary error function is defined by erfc( x ) =√ R∞(2/ π ) x exp(−t2 ) dt. Let rd be a point in the x-y plane with |rd | > 0. We obtainu (rd , t ) =

=

Z t
0

e

Z t
0

Gs (rd , t; 0, s)q(s) ds|rd |2−µ a c(t−s) − 4Dc(t−s)

e

2q(s)

(4πDc(t − s))3/2

"

p
1−

4πDc(t − s) Dc(t2−s)
e `
erfc
2`

p

Dc(t − s)
`

(6)

!#
ds.

This u(rd , t) in Equation (6) is used for the calculation in Section 3.1.
2.1.2. Diffuse Light in Two Dimensions
For the later purpose of a numerical experiment, we here consider light propagation in two
dimensions. Let us suppose that the two-dimensional half space of positive y is occupied witha random medium in the x-y plane. The energy density u of light at position ρ = ( x, y) in the medium
due to the incident beam on the boundary (the x-axis) at ρis (i = 1, . . . , Ms ) obeys the following
diffusion equation.


1 ∂
(7)
− D∆ + µ a (ρ) u(ρ, t; ρis ) = 0,c ∂twith the boundary condition−`∂u(ρ, t; ρis ) + u(ρ, t; ρis ) = δ( x − xsi )δ(t),∂y(8)

Appl. Sci. 2019, 9, 3500

4 of 17and the initial condition u(ρ, 0; ρis ) = 0. The incident laser beam on the right-hand side of theboundary condition was assumed to be a pulse at the position ρis = ( xsi , 0). We suppose that the
diffusion coefficient D is a positive constant but µ a varies in space.
Below we develop the Rytov approximation. Let us write µ a (ρ) ≥ 0 as
µ a (ρ) = µ a0 + δµ a (ρ),(9)

where µ a0 is a constant and the perturbation δµ a (ρ) spatially varies. We note the relation,
Z tZ ∞Z ∞u(ρ, t; ρis ) = u0 (ρ, t; ρis ) −0G (ρ, t; ρ0 , s)δµ a (ρ0 )u(ρ0 , s; ρis ) dx 0 dy0 ds,

−∞

0

(10)

where u0 (ρ, t; ρis ) is the solution to the diffusion Equation (7) with the zeroth-order coefficient replacedby b0 . We introduce the Green’s function G (ρ, t; ρ0 , s), which satisfies

1 ∂− D∆ + µ a0c ∂tG = δ ( ρ − ρ 0 ) δ ( t − s ),(11)with the boundary condition −` ∂G∂y + G = 0 at y = 0, and the initial condition G = 0 at t = 0.
The above relation (10) can be directly verified. By recursive substitution, we obtain u asu(ρ, t; ρis ) = u0 (ρ, t; ρis ) −Z tZ ∞Z ∞00−∞

G (ρ, t; ρ0 , s)δµ a (ρ0 )u0 (ρ0 , s; ρis ) dx 0 dy0 ds + O((δµ a )2 ).

(12)

By neglecting higher-order terms assuming that δµ a is small, we arrive at the (first) Born
approximation [31], in which u is given byu(ρ, t; ρis ) = u0 (ρ, t; ρis ) − c

Z tZ ∞Z ∞0−∞0

G (ρ, t; ρ0 , s)δµ a (ρ0 )u0 (ρ0 , s; ρis ) dx 0 dy0 ds.

(13)

We note that the Green’s function is obtained as [27–30]"( y − y 0 )2
( y + y 0 )2
( x − x 0 )2
e−µa0 c(t−s) − 4Dc−−(t−s) e 4Dc(t−s) + e 4Dc(t−s)
eG (ρ, t; ρ , s) =4πD (t − s)
0

p

−

0 2

(y+y )
4πDc(t − s) − 4Dc(t−s) ee`

y+y0 +2Dc(t−s)/`√
4Dc(t−s)2

erfc

y + y0 + 2Dc(t − s)/`p
4Dc(t − s)

(14)

!#
,

for t > s, and otherwise G (ρ, t; ρ0 , s) = 0. Moreover, we obtaineu0 (ρ, t; ρis ) =−µ a0 ct

2πDt

e

−

( x − xsi )2 +y2
4Dct

√

"
1−



2

`
πDct y+√2Dct/
4Dct
e
erfc
`



y + 2Dct/`
√
4Dct

#
.

(15)

The above expression of u0 is similar to the formula in Equation (6) but does not have a time integral
because in this case we assumed the delta function δ(t) for the incident beam.
To obtain the expression of the Rytov approximation, we introduce ψ0 , ψ1 as [31]u0 = eψ0 ,u = eψ0 +ψ1 .

(16)

Appl. Sci. 2019, 9, 3500

5 of 17

By plugging the above expressions of u, u0 into the Born approximation and neglecting terms ofO((δµ a )2 ), we obtain
Z tZ ∞Z ∞eψ1 = 1 − e−ψ0
G (ρ, t; ρ0 , s)δµ a (ρ0 )u0 (ρ0 , s; ρis ) dx 0 dy0 ds
0 0
−∞


Z tZ ∞Z ∞= exp −e−ψ0
G (ρ, t; ρ0 , s)δµ a (ρ0 )u0 (ρ0 , s; ρis ) dx 0 dy0 ds .
0

(17)

−∞

0

The (first) Rytov approximation is thus derived as"u(ρ, t; ρis ) = u0 (ρ, t; ρis ) exp

1−
u0 (ρ, t; ρis )#Z tZ ∞Z ∞
0

0

−∞

0

0

0

G (ρ, t; ρ , s)δµ a (ρ )u0 (ρ , s; ρis ) dx 0 dy0 ds. (18)Let us define"p( y + y 0 )2
(y+y0 )2 −(y−y0 )24πDc(t − s)
1
−
0
4Dc
(
t
−
s
)
4Dc
(
t
−
s
)
e
g(y, t; y , s) =−1+e
4πD (t − s)
`

2
√
!#
p
0
Dc(t−s)√y+y+Dc(t − s)y + y0`2 Dc(t−s)
p
×e
erfc
+
.
`
2 Dc(t − s)(19)

Then we have
Z t
0

0

0

G (ρ, t; ρ , s)u0 (ρ , s) ds = e−µ a0 ct

Z t
0

0 2

e

0

i 2

(x−x )( x − xs )− 4Dc(t−s) − 4Dcseg(y, t; y0 , s) g(y0 , s; 0, 0) ds.

(20)

Therefore, Equation (18) can be rewritten as"u(ρ, t; ρis ) = u0 (ρ, t; ρis ) exp

×

Z t
0

0 2

e

e−µa0 ct−u0 (ρ, t; ρis )
0

e

0

−∞

δµ a (ρ0 )!

i 2(x−x )( x − xs )− 4Dc(t−s) − 4Dcs

Z ∞Z ∞00g(y, t; y , s) g(y , s; 0, 0) ds

#
0

(21)

0

dx dy .

This expression (21) is used to compute the forward data in Section 3.2.
2.2. Inverse Problems by an Iterative Scheme
We suppose that on the surface of a two- or three-dimensional random medium, for each sourcejjat ρis or ris (i = 1, . . . , Ms ) exiting light is detected at ρd or rd (j = 1, . . . , Md ) and is measured at times
tk (k = 1, . . . , Mt ). Let yijk be measured data whereas u is a solution to the diffusion equation. Let us

suppose that u depends on a vector a which contains unknown parameters. We wish to reconstruct a.
In Section 3.1, a = (µ a , D ), and a is a scalar (a one-dimensional vector) in Section 3.2. For example,
in the former case the solution u depends on a since the calculated value of u depends on µ a , D.
Let us introduce vectors U and F(a) of dimension Ms Md Mt asjjj{F(a)}ijk = u(ρd , tk ; ρis ; a) or u(rd , tk ; ris ; a),{U}ijk = yijk ,
j

j

(22)

j

where we wrote u(ρd , tk ; ρis ) = u(ρd , tk ; ρis ; a) and u(rd , tk ; ris ) = u(rd , tk ; ris ; a) emphasizing that u
depends on a. We find optimal a by minimizing kU − F(a)k22 . Here we particularly consider theLevenberg–Marquardt method [3,4], i.e., the reconstructed value a∗ = limk→∞ ak is computed by theiteration given byhi −1
ak+1 = ak + F 0 (ak ) T F 0 (ak ) + λI
F 0 (ak ) T (U − F(ak )) ,(23)

Appl. Sci. 2019, 9, 3500

6 of 17

where F 0 (a) is the Jacobian matrix, which contains derivatives of F(a) with respect to a, and theparameter λ is nonnegative. By modifying the original algorithm according to Ref. [32], our algorithmof the Levenberg–Marquardt method, which we call Algorithm 1, is described below.
Algorithm 1: Levenberg–Marquardt (LM)
1. Set k = 0 and λ = 1.
2. Calculate F(ak ) and F 0 (ak ).
3. Calculate S(ak ) = dT d, where d = U − F(ak ).
4. Prepare A = F 0 (ak ) T F 0 (ak ) and v = F 0 (ak ) T d.
5. Find δ from ( A + λI )δ = −v.
S(a )−S(a +δ)
6. Obtain S(ak + δ) and R = −δkT (2v+kAδ) .
7. If R < 0.25, then set ν = 10 (αc < 0.1), ν = 1/αc (0.1 ≤ αc ≤ 0.5), or ν = 2 (αc > 0.5), where −1. If R < 0.25 and λ = 0, set λ = 1/k A−1 k and ν = ν/2.
αc = 2 − (S(ak + δ) − S(ak )) /δ T v
In the case of R < 0.25, we set λ = νλ. If R > 0.75, then we set λ = λ/2. If R > 0.75 andλ < 1/k A−1 k, set λ = 0. Otherwise when 0.25 ≤ R ≤ 0.75, no update for λ.
8. If S(ak + δ) ≥ S(ak ), then return to Step 5.
9. If S(ak + δ) < S(ak ), set ak+1 = ak + δ. Then put k + 1 → k and go back to Step 2. Repeat theabove procedure until one of the stopping criteria kδk < 10−4 and S < 10−14 is fulfilled.
2.3. Inverse Problems by Markov Chain Monte Carlo
For simplicity in this section, we describe the Metropolis–Hastings Markov chain Monte Carloalgorithm (MH-MCMC) assuming a is the scalar appearing in Section 3.2. The extension of thederivation to the general case of vector a is straightforward.
Suppose that the coefficient µ a is unknown and µ a depends on a parameter a. We will reconstruct
a within the framework of the Bayesian inversion algorithm [33,34], i.e., we will find the probability
distribution π ( a|U) of a for measured data U. Let f prior ( a) be the prior probability density andπ (U| a) be the likelihood density or the conditional probability density of U for a. Then the joint
probability density π ( a, U) of a, U is given by π ( a, U) = π (U| a) f prior ( a). According to the Bayes
formula, the conditional probability density π ( a|U) is given by
π ( a |U) = R ∞L(U| a) f prior ( a)−∞ L (U| a ) f prior ( a ) da

,

(24)

where L(U| a) is a function proportional to π (U| a). Assuming Gaussian noise, we put
L (U| a ) = e− 12 kU−F( a)k22
2σe

.

(25)

In this paper, we simply set f prior ( a) = 1, i.e., we can say f prior ( a) ∝ 1[ amin ,amax ] ( a) and the interval[ amin , amax ] is large enough so that all a’s appearing in the Markov chain fall into this interval. Here,
1 A ( a) is the indicator function defined as 1 A ( a) = 1 for a ∈ A and = 0 for a ∈/ A. General
uniform distributions can be used for f prior ( a) if we use the prior-reversible proposal that satisfies
f prior ( a)q( a0 | a) = f prior ( a0 )q( a| a0 ) (see below for the proposal distribution q( a0 | a)) [35]. Anotherpossible choice of f prior ( a) is a Gaussian distribution, which turns out to be the Tikhonov regularization
term in the cost function.

Appl. Sci. 2019, 9, 3500

7 of 17

Using the Metropolis–Hastings algorithm, we can evaluate π ( a|U) even when the normalization
factor is not known and only the relation π ( a|U) ∝ L(U| a) f prior ( a) is available [33,34]. We can find
π ( a|U) using a sequence p1 , p2 , . . . as
π ( a|U) = lim pk ( a),

(26)

k→∞

where pk+1 ( a) is obtained from pk ( a) (Markov chain) as
p k +1 ( a 0 ) =Z ∞−∞K ( a0 , a) pk ( a) da.

(27)

For all a, a0 , the transition kernel satisfiesZ ∞K ( a0 , a) ≥ 0,

K ( a0 , a) da0 = 1.

(28)

K ( a 0 , a ) = α ( a 0 , a ) q ( a 0 | a ) + r ( a ) δ ( a 0 − a ),

(29)

−∞

Let us write K ( a0 , a) as
where q( a0 | a) is the proposal distribution andr ( a) = 1 −Z ∞−∞α( a0 , a)q( a0 | a) da0 .

(30)

In the Metropolis–Hastings algorithm we set α( a0 , a) = min {1, π ( a0 |U)q( a| a0 )/[π ( a|U)q( a0 | a)]}.
With this choice of α( a0 , a), the detailed balance is satisfied and K ( a, a0 )π ( a0 |U) = K ( a0 , a)π ( a|U).
A common choice of q( a0 | a) is the normal distribution, i.e., q(·| a) = N ( a, ε2 ) with the mean a and thestandard deviation ε > 0. We note that q( a0 | a) = q( a| a0 ) in this case. We haveZ ∞−∞h( a)π ( a|U) da =

lim

1

kmax

∑ h ( a k ),kmax →∞ k max k =1(31)

where ak ∼ pk (·) and h is an arbitrary continuous bounded function.
Simulated annealing [36] is a type of the Metropolis–Hastings MCMC algorithm in which thetemperature σe decreases during the simulation. The algorithm is summarized below as Algorithm 2.
In this paper we consider two temperatures.
Algorithm 2: Two-temperature simulated annealing (SA)
1. Set large σe as a high temperature.
2. Generate a0 ∼ q(·| ak ) = N ( ak , ε2 ) with ε > 0 for given ak .
3. Calculate α( a0 , ak ) = min {1, π ( a0 |U)/π ( ak |U)}.
4. Update ak as ak+1 = a0 with probability α( a0 , ak ) but otherwise set ak+1 = ak .
5. Continue while k ≤ k b .
6. Then decrease σe to a smaller value as a low temperature, and continue to update ak .
Now we propose an MCMC-iterative hybrid method by combining Algorithms 1 with theMetropolis–Hastings Markov chain Monte Carlo. We first run the Markov chain Monte Carlo.
Although the Metropolis–Hastings MCMC eventually gives the correct solution after about 10,000
steps when the chain reaches the steady state, we stop the Monte Carlo simulation way before theburn-in time. At the k b th Monte Carlo step after initial rapid changes ceases, we record the obtained
reconstructed value akb and switch to Algorithm 1 with the initial value the recorded akb . It is importantthat k b is chosen after the initial rapid changes although k b can be a significant way from the burn-in

Appl. Sci. 2019, 9, 3500

8 of 17

time. Otherwise, the final reconstructed results are quite robust against the choice of k b . We refer to thefollowing hybrid algorithm as Algorithm 3.
Algorithm 3: Hybrid
1. Choose an initial guess a0 , which is not necessarily close to the global minimum.
2. Generate a0 ∼ q(·| ak ) = N ( ak , ε2 ) with ε > 0 for given ak .
3. Calculate α( a0 , ak ) = min {1, π ( a0 |U )/π ( ak |U )}.
4. Update ak as ak+1 = a0 with probability α( a0 , ak ) but otherwise set ak+1 = ak .
5. Obtain a reconstructed akb .
6. Switch to Algorithm 1 with the initial guess akb .
We close this subsection by the Gelman-Rubin convergence diagnostic, which uses the intra-chain
variance and inter-chain variance [37,38]. We run MMC different chains with different initial values.
Let amk denote the kth value in the mth chain (k = 1, . . . , k b , m = 1, . . . , MMC ). We discard the firstk a − 1 steps before the initial rapid changes cease. Then we compute the following intra-chain averageand variance.
kb
kb
1
1
2
ām =
am
,
σ
=
( am − ām )2 .
(32)
∑
∑
m
k b − k a + 1 k=k k
k b − k a k=k kaa

Next we compute the inter-chain mean and variance.
ā =

1 MMC
ām ,
MMC m∑=1

Now we introduceV̂ =where W = M1

B=

k b − k a + 1 MMC( ām − ā)2 .
MMC − 1 m∑
=1kb − k a1W+B,
kb − k a + 1
kb − k a + 1

(33)

(34)

M

2
But W is also an unbiased∑m=MC
1 σm . This is an unbiased estimator of the true variance.p
estimate of the true variance if the chains converge. p
Therefore, we have V̂/W ≈ 1 if converged.
Below, we will see that we can choose k b for which V̂/W is not close to 1. In this paper, we set
MMC = 10.
MC

2.4. TRS Measurements of a Polyurethane-Based Phantom
In this section, we consider time-resolved measurements for a phantom. The solid phantom(biomimic optical phantom) is made of polyurethane to simulate biological tissues (INO, Quebec,
QC, Canada). The absorption coefficient and reduced scattering coefficient are µ a = 0.0209 mm−1
and µ0s = 0.853 mm−1 at wavelength 800 nm. The refractive index of the phantom is n = 1.51.
The measurements were performed by the TRS (time-resolved spectroscopy) instrument (TRS-80,
Hamamatsu Photonics K.K., Hamamatsu, Japan). Two optical fibers from TRS-80 are attached on
the top of the phantom with the separation 13 mm (Ms = Md = 1). The phantom is illuminatedby near-infrared light of the wavelength 760 nm through one optical fiber and the reflected light isdetected by the other optical fiber. The time interval of the measured data is 10 ps and we used the datafrom t1 = 2.00 ns to t Mt = 8.00 ns (Mt = 601). Measured counts, i.e., the number of photons, are shownin Figure 1. The upper panel of Figure 1 shows the instrument response function (IRF), which is given
as a property of the experimental device. The measured reflected light is shown in the lower panel of
Figure 1. The function q(tk ) is the instrument response function divided by the maximum value ofthe measured reflected light. The peak of q is at 2.56 ns. {U }11k = y11k is the measured reflected lightnormalized by its maximum value. In the lower panel of Figure 1, t1 and t Mt are marked by verticaldashed lines. The size of the phantom is large enough that we can assume the three-dimensional halfAppl. Sci. 2019, 9, 3500

9 of 17countsspace. Then the energy density of the detected light is computed by Equation (6). The parameter ` is
obtained from Equation (3).
1.5x1013

IRF

13

1x10

5x1012
0
0

1

2

3

4

5

6

7

8

countstime (ns)1.5x108

reflected light

1x108
5x107
0
0

1

2

3

4

5

6

7

8

time (ns)Figure 1. Measured data from TRS-80. The instrument response function and measured reflected
light are shown in the upper and lower panels, respectively. In the lower panel, the dashed lines show
t1 = 2.00 ns and t Mt = 8.00 ns (Mt = 601).

2.5. Numerical Phantom
To examine when the iterative scheme fails by being trapped by a local minimum and how theMarkov chain Monte Carlo is capable of escaping from such local minima, a toy model is devised
which is simple enough to explicitly understand the structure of the cost function but is complicatedenough that the cost function has one local minimum and one global minimum.
We consider diffuse optical tomography in the two-dimensional space. In our toy model wesuppose that the diffusion coefficient D is constant everywhere in the medium but there is absorption
inhomogeneity and the absorption coefficient δµ a is unknown. Moreover, we assume that δµ a (ρ) is
given byδµ a (ρ) = η f a ( x )δ(y − y0 ),(35)
where η, y0 are given positive constants. Here, f a ( x ) is given bytanh x2
f a (x) = a + 3 1 +
10


3




a

2




1 − tanh x2 ,(36)

where a > 1.1 is a constant. Thus, a is the parameter to be reconstructed in our toy inverse problem.
As is shown in Figure 2, the function f a ( x ) has one peak at x = 0 and the maximum value isf a (0) = a2 ( a + 3), f a monotonically decays for | x | > 0, and f a → 0 as | x | → ∞.
Now we describe how the forward data is computed. The unit of length and unit of time aretaken to be mm and ps, respectively. On the x-axis, we place two sources (Ms = 2) at ρ1s = (−20, 0),ρ2s = (20, 0) and three detectors (Md = 3) at ρ1d = (−40, 0), ρ2d = (0, 0), ρ3d = (40, 0). We set µ0s = 1,
µ a = 0.02, n = 1.37. Suppose that there is absorption inhomogeneity at depth 5. For δµ a , we putη = 300/c, y0 = 5, anda = 1.5.
(37)
To distinguish, hereafter the true value of a is denoted by ā. We assume 3% Gaussian noise and givethe measured data asjyijk = u(ρd , tk ; ρis ; ā)(1 + eijk ),(38)
where eijk ∼ N (0, 0.032 ). In this numerical experiment, Mt = 500 and tk+1 − tk = t1 = 5(k = 1, . . . , Mt − 1).

Appl. Sci. 2019, 9, 3500

10 of 17

16

a=1.7
a=1.5
a=1.3

14
12

fa

10
8
6
4
2
0
-4

-2

0

2

4

x
Figure 2. The function f a ( x ) in Equation (36) is plotted for a = 1.3, 1.5, and 1.7.

By substituting the form (35) of δb in (21), we can express the energy density asjju(ρd , tk ; ρis ; a) = u(ρd , tk ; ρis )"j= u0 (ρd , tk ; ρis ) exp

× g(y0 , s; 0, 0) 
where
u0 =

k
e−µa0 ct

2πDtk

Z ∞
−∞

√



j
( x − x i )2
− d sk 1
e 4Dct−−ηe−µa0 ct

Z tkju0 (ρd , tk ; ρis )f a ( x 0 )e

πDctk

`

k

j
( x − x 0 )2
− d k
4Dc(t −s)

√

e

Dctk
`

0

g(0, tk ; y0 , s)e−

( x 0 − xsi )2
4Dcs
dx 0  ds

√

2

erfc

Dctk

`

!
.

(39)
#
,

(40)

3. Results
3.1. Determination of Optical Properties
Here we consider reconstructions from measured data in the phantom experiment. Let us firstconsider initial guesses below.
Case 1:µ a = 0.01 mm−1 ,µ0s = 1.0 mm−1 .

(41)

In this Case 1, the following µ a , µ0s are obtained both by Algorithm 1 (LM) and Algorithm 3 (hybrid).
µ a = 0.016 mm−1 ,µ0s = 0.63 mm−1 .

(42)

The results for µ a and µ0s are shown in Figures 3 and 4, respectively. In Figure 3, reconstructed valuesof µ a are plotted against the iteration number k. In the top panel of Figure 3, we see that Algorithm 1(LM) quickly converges. In Algorithms 2 and 3, we set k b = 99. As is seen in the middle panel of
Figure 3, the convergence of Algorithm 2 (SA) is slow. The temperature is decreased from σe = 10−6 toσe = 10−7 at the k b th Monte Carlo step. Similarly, ε is changed
p from 0.1 to 0.001. Algorithm 2 returns
the correct values after many Monte Carlo steps; we have V̂/W = 1.02, 1.06 for µ a , µ0s , respectively,
when k a = 9000 and k b = 10000. Finally, in the bottom panel of Figure 3, we see that the iteration of
Algorithm 3 (hybrid) immediately converges after we switch from the Metropolis–Hastings MCMC
−6(σe = 10p , ε = 0.1) to Algorithm 1 (LM). No convergence of the Monte Carlo chain is required, and we
found V̂/W = 8.6 for µ a and = 11.0 for µ0s (k a = 49, k b = 99). A similar behavior is observed for thereconstruction of µ0s in Figure 4. In the top panel of Figure 4, we find that Algorithm 1 (LM) works best

Appl. Sci. 2019, 9, 3500

11 of 17and converges after a few iterations, whereas Algorithm 2 (SA) in the middle panel of Figure 4 has slow
convergence and reconstructed values around at k b = 99 are still away from µ0s in (42). By switchingfrom MH-MCMC to Algorithm 1 (LM) using Algorithm 3 (hybrid), convergence is easily obtained asshown in the bottom panel of Figure 4.
Now we start the simulation by setting the following initial values.
µ a = 0.5 mm−1 ,Case 2:µ0s = 1.0 mm−1 .

(43)µa (mm-1)

In this Case 2, Algorithm 1 (LM) fails and returns µ a = 0.068 mm−1 and µ0s = 1.75 mm−1 whereas
Algorithm 3 (hybrid) still gives the correct values. The reconstructed values of µ a and µ0s at each
iteration are shown in Figures 5 and 6, respectively. In Figure 5, the vertical axes of the left three panels
are from 0 to 0.6 whereas the vertical axes of the right three panels are from 0 to 0.08. In Figure 6,
the vertical axes of the left three panels are from 0 to 2 while the vertical axes of the right three panels
are from 0 to 1.3. In the top panels of Figure 5, we see that the reconstruction of µ a is unsuccessful
by Algorithm 1 (LM). Algorithm 2 (SA) approaches µ a in Equation (42) but still deviates from thatvalue in the middle panels of Figure 5. As shown in the bottom panels of Figure 5, Algorithm 3(hybrid) converges in a few iterations after switching to Algorithm 1 (LM). In the top panels of Figure 6,
Algorithm 1 (LM) fails to reconstruct µ0s . In the middle panels of Figure 6, reconstructed values by
Algorithm 2 (SA) come close to µ0s in Equation (42) but suffer from slow convergence. In the bottom
panels of Figure 6, we see that Algorithm 3 (hybrid) arrives at µ0s in Equation (42) without any problem.
0.08
0.06
0.04
0.02
0

Levenberg-Marquardt
0

20

40

60

80

100

120

µa (mm-1)

k
0.08
0.06
0.04
0.02
0

simulated annealing
0

20

40

60

80

100

120

80

100

120

µa (mm-1)

k
0.08
0.06
0.04
0.02
0

hybrid
0

20

40

60

k

µs’ (mm-1)

Figure 3. Case 1: Reconstructed µ a by (top) Algorithm 1, (middle) Algorithm 2, and (bottom)
Algorithm 3. The dashed lines show µ a in Equation (42).
1.2
0.8
0.4
0

Levenberg-Marquardt
0

20

40

60

80

100

120

µs’ (mm-1)

k
1.2
0.8
0.4
0

simulated annealing
0

20

40

60

80

100

120

80

100

120

µs’ (mm-1)

k
1.2
0.8
0.4
0

hybrid
0

20

40

60

k

Figure 4. Case 1: Reconstructed µ0s by (top) Algorithm 1, (middle) Algorithm 2, and (bottom)
Algorithm 3. The dashed lines show µ0s in Equation (42).

12 of 17

0.6
0.4µa (mm-1)µa (mm-1)

Appl. Sci. 2019, 9, 3500

Levenberg-Marquardt

0.2
0
0

20

40

60

80

100

0.08
0.06
0.04
0.02
0

120

Levenberg-Marquardt
0

20

40

60

0.6
0.4

simulated annealing

0.2
0
0

20

40

60

80

100

0.08
0.06
0.04
0.02
0

120

0

µa (mm-1)µa (mm-1)

hybrid

0.2
0
20

40

60

120

20

40

60

80

100

120

80

100

120

k

0.6

0

100

simulated annealing

k
0.4

80

k
µa (mm-1)µa (mm-1)

k

80

100

0.08
0.06
0.04
0.02
0

120

hybrid
0

20

40

60

k

k

2
1.6
1.2
0.8
0.4
0

µs’ (mm-1)µs’ (mm-1)

Figure 5. (Left) Case 2: Reconstructed µ a by (top) Algorithm 1, (middle) Algorithm 2, and (bottom)
Algorithm 3. The dashed lines show µ a in Equation (42). (Right) Same as the left panel but the verticalaxes are from 0 to 0.08.
Levenberg-Marquardt

0

20

40

60

80

100

1.2
0.8
0.4
0

120

Levenberg-Marquardt
0

20

40

60

2
1.6
1.2
0.8
0.4
0

simulated annealing

0

20

40

60

80

100

1.2
0.8
0.4
0

120

20

40

60

0

µs’ (mm-1)µs’ (mm-1)
80

k

120

20

40

60

80

100

120

80

100

120

k

hybrid0100simulated annealing

k
2
1.6
1.2
0.8
0.4
0

80

k
µs’ (mm-1)µs’ (mm-1)

k

100

120

1.2
0.8
0.4
0

hybrid
0

20

40

60

k

Figure 6. (Left) Case 2: Reconstructed µ0s by (top) Algorithm 1, (middle) Algorithm 2, and (bottom)
Algorithm 3. The dashed lines show µ0s in Equation (42). (Right) Same as the left panel but the verticalaxes are from 0 to 1.3.

The initial guesses for Case 1 are reasonably close to the values found in Equation (42). However,
we started with initial guesses which are quite different from the above-mentioned values in Case 2.
It is not surprising that Algorithm 1 (LM) does not work for Case 2, but Algorithm 3 (hybrid) canarrive at the correct values. The numerical calculations were performed on Matlab (i7-8700 CPU, 16 GB
memory). In the hybrid method, the Metropolis–Hastings MCMC does not reach the steady stateat about 100 steps but can be switched to the Levenberg–Marquardt method. For Figures 5 and 6,
the computation time was 5 s. The simulated annealing method returns the correct value after about
10,000 steps, but it takes 8 min. The computation for Algorithm 1 (LM) stopped in 0.5 s.
Below we summarize the reconstructed values on Table 1. Although Algorithm 2 (SA) and
Algorithm 3 (hybrid) return the same results after a long time, the hybrid scheme is about ten times
faster, and moreover there is no need for choosing the low temperature for the latter algorithm.
For Algorithm 3, it is not necessary to wait until the burn-in time, but it is enough if the initial rapidchange ceases. In Section 3.2, we illustrate that Algorithm 3 works once the Monte Carlo simulation
escapes from a local minimum and the algorithm does not require that the calculation reaches thesteady state.

Appl. Sci. 2019, 9, 3500

13 of 17

Table 1. Reconstructed values of µ a , µ0s are shown for Case 1 and Case 2. The units of µ a , µ0s are bothmm−1 . For Algorithm 2, values at 120th Monte Carlo step are shown.

initial values
Algorithm 1 (LM)
Algorithm 2 (SA)Algorithm 3 (hybrid)

Case 1 (µ a , µ0s )Case 2 (µ a , µ0s )(0.01, 1.0)
(0.016, 0.63)(0.032, 1.02)(0.016, 0.63)(0.5, 1.0)
(0.068, 1.75)(0.028, 0.92)(0.016, 0.63)

3.2. Determination of Absorption Inhomogeneity
We perform a numerical experiment of diffuse optical tomography using the toy model. Let us∂
consider when the iterative scheme fails. We see F 0 (0) = ∂au ( ρ j , t k ; ρi ; a )= 0. For sufficiently

a =0
small e > 0, we have F 0 (e) > 0, F 0 (−e) < 0, | F 0 (±e)| ' 0, and U − F (±e) < 0. Therefore, if we start

the iteration from a0 = e, we obtain| a1 | < e,

| a2 | < | a1 |,

....

(44)

Thus, the sequence ak approaches 0 and can never arrive at ā (= 1.5).
Let us consider how the difference U − F ( a) depends on a. We introducey2
1h(t) = exp − 0
t
4Dct

√

!"
1−

√

0 + Dct
πDct 2√yDct
`
e
`



2

erfc

y
√0 +
2 Dct

√

Dct
`

!#
.

(45)

The following form is obtained using Equation (39). By neglecting noise, we haveu(ρ j , tk ; ρi ; ā) − u(ρ j , tk ; ρi ; a)k Z

ηe−µa0 ct=(2πD )2tk
0
h(tk − s) h(s) 

Z ∞
−∞


j
 − ( x d − x 0 )2
0 − x i )2(xs2d ā ( a, x 0 ) 1 − tanh x 0 e 4Dc(tk −s) e− 4Dcs dx 0  ds,

(46)

where d ā ( a, x 0 ) = ξ ( ā, x 0 ) − ξ ( a, x 0 ) with"0

ξ ( a, x ) = a

2tanh x 0
a+3 1+
10

2

!#
.

(47)

0
For a given|d ā ( a; x 0 )|2 has one global minimum at a = ā, one local minimum at
 x , the function02x
a = −2 1 + tanh
, and one local maximum at a = 0. Therefore, the above expression implies
10
that Algorithm 1 (LM) fails when the initial value a0 is negative and the correct value ā = 1.5 isreconstructed only for a positive initial guess. Indeed, the computation ends up with −2.05 if westart the iteration from a0 = −0.1 (see below) or −0.01, and the value 1.68 is obtained when a0 = 0.01.
Figure 7 shows |d ā ( a, x )|2 for ā = 1.5 and tanh( x2 ) = 0.5.
In Figure 8, we plot reconstructed values of a against iteration numbers. The initial value is set
to a0 = −0.1. The reconstruction by Algorithm 1 (LM) is shown in the top panel of Figure 8. As we
can predict from Figure 7, Algorithm 1 (LM) converges to the local minimum and can never arriveat the global minimum. Monte Carlo simulation can jump over the local maximum and approach
the global minimum, but keeps fluctuating as shown in the middle panel of Figure 8. We initiallyset (σe , ε) = (10−6 , 0.5) for Algorithms 2 and 3. In Algorithm 2, we set (σe , ε) = (10−7 , 0.005) after
the temperature decreases at the k b th Monte Carlo step. In the bottom panel of Figure 8, Algorithm 3(hybrid) successfully arrives at the global minimum. Using Matlab, the computation time was 17 min
while we needed 3 hr for Algorithm 2 (SA) (1000 steps).

Appl. Sci. 2019, 9, 3500

14 of 17160

2|da- (a,x)|
a- =1.5 and tanh(x2)=0.5

140
120
100
80
60
40
20
0
-3

-2

-1

a

0

1

2

Figure 7. The function |d ā ( a, x )|2 (ā = 1.5) is plotted when tanh( x2 ) = 0.5.

aAfter the initial time with large jumps, it is found that we can set k b = 99 in our simulation of the
MCMC-iterative hybrid method. The correct value of ā is reconstructed by Algorithm 3 (hybrid) evenwhen the simulation experiences a local minimum. Starting from akb = 1.8257, the calculation stops ata∗ = 1.6841 (a∗ = akb +3 ). We note that the reconstructed a∗ is not exactly 1.5 due to noise.
When the initial guess a0 = −0.1 lies in the valley of the local minimum (see Figure 7),
the reconstructed a by Algorithm 3 (hybrid) moves to the valley of the global minimum with thehelp of Monte Carlo simulation as shown in Figure 8, while the reconstructed a by Algorithm 1 (LM)falls to the local minimum as the nature of Newton-type iterative methods. For Algorithm 3 (hybrid),
we obtained a∗ = 1.6841. There is a possibility that negative reconstructed values are obtained by
Algorithm 2 (SA) and Algorithm 3 (hybrid) for the first a few iterations if different pseudo-random
numbers are used. These negative reconstructed values, however, will turn positive and the behaviorof reconstructed a by Algorithm 2 (SA) and Algorithm 3 (hybrid) is always more or less similar to the
middle and bottom panels of Figure 8. For the initial guesses of a0 = −0.01 and a0 = 0.01, Algorithm 3(hybrid) works without any problem and Algorithm 2 (SA) also returns a reasonable result after
a sufficiently large number of iterations (results not shown).
2
1
0
-1
-2

Levenberg-Marquardt
0

20

40

60

80

100

120

80

100

120

80

100

120

a

k
2
1
0
-1
-2

simulated annealing
0

20

40

60

a

k
2
1
0
-1
-2

hybrid
0

20

40

60

k
Figure 8. Reconstructed ā (= 1.5) by (top) Algorithm 1 (LM), (middle) Algorithm 2 (SA), and (bottom)
Algorithm 3 (hybrid) for the initial value a0 = −0.1.

4. Discussion
In this paper, we have proposed a hybrid numerical scheme which uses Markov chain Monte Carloin the first step and then uses an iterative method in the second step. We switch from MH-MCMC to
LM by observing proposed parameter values. For the typical jump length of parameters in MH-MCMC,

Appl. Sci. 2019, 9, 3500

15 of 17ε = 0.1 was used in Section 3.1 and ε = 0.5 was chosen in Section 3.2. Although these values were
set so that the MH-MCMC calculation was efficiently performed, other choices of ε are also possible.
The proposed scheme is quite general and can be applied to different inverse problems in NIRS which
are solved by iterative methods even when the forward problem must be solved fully numericallywith finite difference method or finite element method. It is an interesting future study how the hybrid
scheme can be extended to diffuse optical tomography, which has many unknowns.
More sophisticated algorithms than the Metropolis–Hastings Markov chain Monte Carlo usedin this paper have been proposed to overcome slow convergence. The delayed rejection scheme
reduces the net rejection rate [39]. In the adaptive Metropolis algorithm, parameters in the proposal
distribution are adjusted during Monte Carlo steps [40]. The DRAM algorithm, which combines theabove-mentioned two schemes, was also proposed [41]. Two-level MCMC algorithms [42,43] anda multi-level MCMC [23] have been investigated to improve the MCMC algorithm. Such Monte
Carlo schemes might improve the first step of our hybrid method by finding the valley of the global
minimum more easily.
Related to the Metropolis–Hastings Markov chain Monte Carlo algorithm, quantum annealing [44]has been developed in addition to simulated annealing. Aiming at escaping from local minima,brute-force search and genetic algorithm [45] are also well-known optimization algorithms.
The introduction of these methods in NIRS may be found useful in the future.
For the clinical use of NIRS, it has been recognized from early days that finding absolute valuesof the absorption and scattering coefficients is important [2]. In Ref. [15], it is emphasized that theobtained absolute values highly depend on the starting values of the initial estimate for the study ofthe infant brain with their measurement system. By performing Markov chain Monte Carlo before
using standard iterative schemes, we may automatically acquire good initial values. Such clinicalapplication is a natural next step of our research.
Author Contributions: Conceptualization, M.M. and G.N.; methodology, Y.J. and M.M.; software, Y.J. and M.M.;
formal analysis, Y.J.; investigation, Y.H.; data curation, Y.J. and Y.H.; writing—original draft preparation, M.M.
Funding: The first author (Y.J.) is supported by National Natural Science Foundation of China (No. 11971121).
Y.H. and M.M. acknowledge support from Grant-in-Aid for Scientific Research (17H02081) of the Japan Society for
the Promotion of Science. M.M. also acknowledges support from Grant-in-Aid for Scientific Research (17K05572)
of the Japan Society for the Promotion of Science and from the JSPS A3 foresight program: Modeling andComputation of Applied Inverse Problems. G.N. was supported by Grant-in-Aid for Scientific Research (15K21766
and 15H05740) of the Japan Society for the Promotion of Science.
Acknowledgments: The authors appreciate Goro Nishimura for fruitful discussion on optical tomography
and the use of his computer facilities at Hokkaido University. We thank Tatsunori Emi, Takeshi Iwasaki,
and Tomoya Sugiyama for helping related experiments which stimulated the present work.
Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of thestudy; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to
publish the results.

References
1.
2.
3.
4.
5.

Boas, D.A.; Brooks, D.H.; Miller, E.L.; DiMarzio, C.A.; Kilmer, M.; Gaudette, R.J.; Zhang, Q. Imaging the
body with diffuse optical tomography. IEEE Signal Process. Mag. 2001, 18, 57–75. [CrossRef]
Hoshi, Y. Towards the next generation of near-infrared spectroscopy. Phil. Trans. R. Soc. A 2011, 369,
4425–4439. [CrossRef] [PubMed]
Levenberg, K. A method for the solution of certain non-linear problems in least squares. Q. Appl. Math. 1944,
2, 164–168. [CrossRef]
Marquardt, D.W. An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math.
1963, 11, 431–441. [CrossRef]
Lange, F.; Tachtsidis, I. Clinical brain monitoring with time domain NIRS: A review and future perspectives.
Appl. Sci. 2019, 9, 1612. [CrossRef]

Appl. Sci. 2019, 9, 3500

6.
7.
8.
9.

10.

11.
12.
13.
14.

15.

16.
17.
18.

19.
20.
21.

22.
23.

24.

25.
26.
27.
28.
29.

16 of 17

Chance, B.; Zhuang, Z.; Unah, C.; Alter, C.; Lipton, L. Cognition-activated low frequency modulation of lightabsorption in human brain. Proc. Natl. Acad. Sci. USA 1993, 90, 3770–3774. [CrossRef]
Hoshi, Y.; Tamura, M. Detection of dynamic changes in cerebral oxygenation coupled to neuronal function
during mental work in man. Neurosci. Lett. 1993, 150, 5–8. [CrossRef]
Kato, T.; Kamei, A.; Takashima, S.; Ozaki, T. Human visual cortical function during photic stimulate on
monitoring by means of near-infrared spectroscopy. J. Cereb. Blood Flow Metab. 1993, 13, 516–520. [CrossRef]
Villringer, A.; Plank, J.; Hock, C.; Schleikofer, L.; Dirnagl, U. Near-infrared spectroscopy (NIRS): A new tool
to study hemodynamic changes during activation of brain function in human adults. Neurosci. Lett. 1993,
154, 101–104. [CrossRef]
Laidevant, A.; da Silva, A.; Berger, M.; Dinten, J.-M. Effects of the surface boundary on the determinationof the optical properties of a turbid medium withtime-resolved reflectance. Appl. Opt. 2006, 45, 4756–4764.
[CrossRef]
Gibson, A.P.; Hebden, J.C.; Arridge, S.R. Recent advances in diffuse optical imaging. Phys. Med. Biol. 2005,
50, R1–R43. [CrossRef]
Arridge, S.R. Methods in diffuse optical imaging. Philos. Trans. R. Soc. A 2011, 369, 4558–4576. [CrossRef]
Arridge, S.R.; Schweiger, M. A gradient-based optimization scheme for optical tomography. Opt. Express
1998, 2, 213–226. [CrossRef]
Choe, R.; Corlu, A.; Lee, K.; Durduran, T.; Konecky, S.D.; Grosicka-Koptyra, M.; Arridge, S.R.; Czerniecki,
B.J.; Fraker, D.L.; DeMichele, A.; et al. Diffuse optical tomography of breast cancer during neoadjuvant
chemotherapy: A case study with comparison to MRI. Med. Phys. 2005, 32, 1128–1139. [CrossRef]
Hebden, J.C.; Gibson, A.; Austin, T.; Yusof, R.M.; Everdell, N.; Delpy, D.T.; Arridge, S.R.; Meek, J.H.; Wyatt,
J.S. Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional
optical tomography. Phys. Med. Biol. 2004, 49, 1117–1130. [CrossRef]
Arridge, S.R.; Schweiger, M. Image reconstruction in optical tomography. Philos. Trans. R. Soc. Lond. B 1997,
352, 717–726. [CrossRef]
Schweiger, M.; Arridge, S. The toast++ software suite for forward and inverse modeling in optical
tomography. J. Biomed. Opt. 2014, 19, 040801. [CrossRef]
Zhao, H.; Gao, F.; Tanikawa, Y.; Homma, K.; Yamada, Y. Time-resolved diffuse optical tomographic imaging
for the provision of both anatomical and functional information about biological tissue. Appl. Opt. 2005, 44,
1905–1916. [CrossRef]
Arridge, S.R. Optical tomography in medical imaging. Inverse Probl. 1999, 15, R41–R93. [CrossRef]
Arridge, S.R.; Schotland, J.C. Optical tomography: Forward and inverse problems. Inverse Probl. 2009, 25,
123010. [CrossRef]
Arridge, S.R.; Kaipio, J.P.; Kolehmainen, V.; Schweiger, M.; Somersalo, E.; Tarvainen, T.; Vauhkonen, M.
Approximation errors and model reduction with an application in optical diffusion tomography. Inverse Probl.
2006, 22, 175–195. [CrossRef]
Bal, G.; Langmore, I.; Marzouk, Y. Bayesian inverse problems with monte carlo forward models.
Inverse Probl. Imaging 2013, 7, 81–105.
Langmore, I.; Davis, A.B.; Bal, G. Multipixel retrieval of structural and optical parameters in a 2-d scene
with a path-recycling monte carlo forward model and a new bayesian inference engine. IEEE Trans. Geosci.
Remote Sens. 2013, 51, 2903–2919. [CrossRef]
Bamett, A.H.; Culver, J.P.; Sorensen, A.G.; Dale, A.; Boas, D.A. Robust inference of baseline optical propertiesof the human head with three-dimensional segmentation from magnetic resonance imaging. Appl. Opt. 2003
42, 3095–3108.
Groenhuis, R.A.J.; Ferwerda, H.A.; Ten Bosch, J.J. Scattering and absorption of turbid materials determinedfrom reflection measurements. 1: Theory. Appl. Opt. 1983, 22, 2456–2462. [CrossRef]
Egan, W.G.; Hilgeman, T.W. Optical Properties of Inhomogeneous Materials; Academic Press: New York, NY,
USA, 1979.
Carslaw, H.S.; Jaeger, J.C. Conduction of Heat in Solids; Oxford University Press: London, UK, 1959.
Yosida, K.; Ito, S. Functional Analysis and Differential Equations; Iwanami: Tokyo, Japan, 1976. (In Japanese)
Hielscher, A.H.; Jacques, S.L.; Wang, L.; Tittel, F.K. The influence of boundary conditions on the accuracy ofdiffusion theory in time-resolved reflectance spectroscopy of biological tissues. Phys. Med. Biol. 1995, 40,
1957–1975. [CrossRef]

Appl. Sci. 2019, 9, 3500

30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.

17 of 17

Machida, M.; Nakamura, G. Born series for the photon diffusion equation perturbing the Robin boundary
condition. arXiv 2017, arXiv:1706.04500.
Ishimaru, A. Wave Propagation and Scattering in Random Media; Academic Press: New York, NY, USA, 1978.
Fletcher, R. A Modified Marquardt Subroutine for Nonlinear Least Squares (Report AERE-R 6799); The Atomic
Energy Research Establishment: Harwell, UK, 1971.
Kaipio, J.; Somersalo, E. Statistical and Computational Inverse Problems; Springer: New York, NY, USA, 2005.
Nakamura, G.; Potthast, R. Inverse Modeling; IOP Publishing: Bristol, UK, 2015.
Iglesias, M.A.; Lin, K.; Stuart, A.M. Well-posed bayesian geometric inverse problems arising in subsurface
flow. Inverse Probl. 2014, 30, 114001. [CrossRef]
Kirkpatrick, S.; Gelatt, C.D., Jr.; Vecchi, M.P. Optimization by simulated annealing. Science 1983, 220, 671–680.
[CrossRef]
Gelman, A.; Rubin, D.B. Inference from iterative simulation using multiple sequences. Stat. Sci. 1992 7,
457–472. [CrossRef]
Martinez W.L.; Martinez, A.R. Computational Statistics Handbook with MATLAB; Chapman and Hall/CRC:
London, UK, 2015.
Tierney, L. Markov chains for exploring posterior distributions. Ann. Stat. 1994, 22, 1701–1762. [CrossRef]
Haario, H.; Saksman, E.; Tamminen, J. An adaptive metropolis algorithm. Bernoulli 2001, 7, 223–242.
[CrossRef]
Haario, H.; Laine, M.; Mira, A.; Saksman, E. Dram: Efficient adaptive mcmc. Stat. Comput. 2006, 16, 339–354.
[CrossRef]
Christen, J.A.; Fox, C. Markov chain monte carlo using an approximation. J. Comput. Graph. Stat. 2005, 14,795–810. [CrossRef]
Efendiev, Y.; Hou, T.; Luo, W. Preconditioning markov chain monte carlo simulations using coarse-scalemodels. SIAM J. Sci. Comput. 2006, 28, 776–803. [CrossRef]
Kadowaki, T.; Nishimori, H. Quantum annealing in the transverse ising model. Phys. Rev. E 1998, 58,
5355–5363. [CrossRef]
Holland, J.H. Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology,
Control and Artificial Intelligence; University of Michigan Press: Ann Arbor, MI, USA, 1975.
c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open accessarticle distributed under the terms and conditions of the Creative Commons Attribution(CC BY) license (http://creativecommons.org/licenses/by/4.0/).