Abstract

Diffuse optical tomography (DOT) is an imaging modality that uses near-infrared light. Although iterative numerical schemes are commonly used for its inverse problem, correct solutions are not obtained unless good initial guesses are chosen. We propose a numerical scheme of DOT, which works even when good initial guesses of optical parameters are not available. We use simulated annealing (SA), which is a method of the Markov-chain Monte Carlo. To implement SA for DOT, a spin Hamiltonian is introduced in the cost function, and the Metropolis algorithm or single-component Metropolis–Hastings algorithm is used. By numerical experiments, it is shown that an initial random spin configuration is brought to a converged configuration by SA, and targets in the medium are reconstructed. The proposed numerical method solves the inverse problem for DOT by finding the ground state of a spin Hamiltonian with SA.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. INTRODUCTION

In near-infrared spectroscopy, tomographic images of optical properties are obtained by diffuse optical tomography (DOT) [1,2]. To obtain reconstructed images, inverse problems of determining coefficients of the diffusion equation from boundary measurements are solved.

For these inverse problems in DOT, direct approaches such as the use of the Born approximation and iterative methods such as the conjugate gradient method and damped Gauss–Newton method are commonly used [3,4]. One of the direct approaches is the inversion of a linearized problem by singular value decomposition (SVD) with the ${L^2}$ regularization (for example, see [5]). To solve the nonlinear inverse problem of DOT, iterative methods are frequently used [3]. For example, breast cancer was detected with the nonlinear-conjugate gradient method [6]. See, e.g., [7,8] for recent advances in DOT.

Although different penalty terms can be employed in iterative methods, the calculation gets trapped by local minima unless a good initial guess is given (for example, see [9]). It was demonstrated in [10] that bad initial guesses (even though they look reasonable) lead to completely wrong reconstructed results, even when the diffusion coefficient is known, and the absorption coefficient is characterized by only one unknown parameter in the diffusion equation.

As numerical schemes that are free from the issue of being trapped by local minima, statistical approaches have been developed [11,12]. Monte Carlo methods have the advantage that the computation is not trapped by local minima, since jumps occur randomly in the landscape. In particular, the Metropolis–Hastings Monte Carlo algorithm was used to estimate parameters in coefficients of the diffusion equation or the radiative transport equation, which is approximated to the diffusion equation at large scales [10,13,14]. However, only several unknown parameters, such as the position of a target and the absorption and scattering coefficients of the target, were determined in those studies.

In this paper, values of the absorption coefficient will be determined at 1830 different points in the medium. The simulated annealing (SA) (e.g., [15]) makes it possible to handle a large number of unknown parameters. To the best of our knowledge, this is the first attempt to use SA for DOT.

The rest of the paper is organized as follows. In Section 2, we introduce the diffusion equation, in which near-infrared light in biological tissue is obeyed. In Section 3, we formulate our inverse problem with the Rytov approximation. Spins are introduced in Section 4. Then, the cost function for the inverse problem is rewritten as a spin Hamiltonian in Sections 5 and 6. The algorithm of SA for DOT is described in Section 7. In Section 8, numerical tests of our SA are presented. Finally, concluding remarks are given in Section 9.

2. DIFFUSION EQUATION

We consider DOT in a domain $\Omega \subset {{\mathbb{R}}^{d}}$ ($d = 2,3$). Let $\partial \Omega$ be the boundary of $\Omega$. Let ${\boldsymbol \nu}({\boldsymbol r})$ be the outer unit normal vector at ${\boldsymbol r} \in \partial \Omega$. We assume that $\Omega$ is occupied by biological tissue, and the outside of $\Omega$ is vacuum. Let $u$ be the photon density of near-infrared light.

If the reconstruction is done using time-resolved data, the time-dependent diffusion equation needs to be considered. Let $c$ be the speed of light in $\Omega$. We assume that the diffusion coefficient ${D_0}$ is a positive constant but ${\mu _a}$ varies in space. The energy density $u$ obeys

$$\left\{{\begin{array}{*{20}{l}}{\left({\frac{1}{c}\frac{\partial}{{\partial t}} - {D_0}\Delta + {\mu _a}({\boldsymbol r})} \right)u({\boldsymbol r},t) = f({\boldsymbol r},t),\quad {\boldsymbol r} \in \Omega ,\quad t \gt 0,}\\[4pt]{{D_0}{\boldsymbol \nu} \cdot \nabla u({\boldsymbol r},t) + \frac{1}{\zeta}u({\boldsymbol r},t) = 0,\quad {\boldsymbol r} \in \partial \Omega ,\quad t \gt 0,}\\[4pt]{u({\boldsymbol r},0) = 0,\quad {\boldsymbol r} \in \Omega ,}\end{array}} \right.$$
where $\zeta \gt 0$ is a constant (see Section 8), and $f({\boldsymbol r},t)$ is the source term. We note that ${\mu _a}$ is non-negative.

In the case of the continuous-wave measurement, light propagation in $\Omega$ is governed by the following diffusion equation for $u({\boldsymbol r})$ (${\boldsymbol r} \in \Omega$):

$$\left\{{\begin{array}{*{20}{l}}{- {D_0}\Delta u({\boldsymbol r}) + {\mu _a}({\boldsymbol r})u({\boldsymbol r}) = f({\boldsymbol r}),\quad {\boldsymbol r} \in \Omega ,}\\[4pt]{{D_0}{\boldsymbol \nu} \cdot \nabla u({\boldsymbol r}) + \frac{1}{\zeta}u({\boldsymbol r}) = 0,\quad {\boldsymbol r} \in \partial \Omega .}\end{array}} \right.$$

Even for measurements in time domain, quite often the time-independent diffusion of Eq. (2) is used for reconstruction with moments or Laplace transforms of the raw time-resolved data. We assume the incident beam $f({\boldsymbol r})$ as

$$f({\boldsymbol r},t) = {g_0}\delta ({\boldsymbol r} - {\boldsymbol r}_s^{(p)})\delta (t)\quad {\rm or}\quad f({\boldsymbol r}) = {g_0}\delta ({\boldsymbol r} - {\boldsymbol r}_s^{(p)}),$$
where ${g_0} \gt 0$ is a constant, and ${\boldsymbol r}_s^{(p)}$ is the position of the source of the $p$th source-detector pair ($p = 1,2, \ldots ,{M_{{\rm SD}}}$). Furthermore, $\delta (\cdot)$ is the Dirac delta function. Light is detected on the boundary at ${\boldsymbol r}_d^{(p)}$.

We choose a region of interest ${\Omega _{{\rm ROI}}}$ in $\Omega$. Let ${\bar \mu_a}$ be a non-negative constant. We suppose

$${\mu _a}({\boldsymbol r}) = {\bar \mu_a}$$
for ${\boldsymbol r}$ in $\Omega$, including the boundary $\partial \Omega$ and excluding ${\Omega _{{\rm ROI}}}$. Let us write ${\mu _a}({\boldsymbol r})$ as
$${\mu _a}({\boldsymbol r}) = {\bar\mu_a} + \delta {\mu _a}({\boldsymbol r}),\quad {\boldsymbol r} \in {\Omega _{{\rm ROI}}}.$$

3. RYTOV APPROXIMATION

We will use SA by regarding $\delta {\mu _a}({\boldsymbol r})$ as spins. For this purpose, here, we develop the perturbation theory and give the solution $u({\boldsymbol r})$ in terms of an integral of $\delta {\mu _a}({\boldsymbol r})$. Since the Rytov approximation for the time-dependent case is similarly derived, we consider the time-independent case below.

Let ${u_0}({\boldsymbol r})$ be the solution to Eq. (2) in the case of $\delta {\mu _a} \equiv 0$. Then, ${u_0}({\boldsymbol r}) = {g_0}G({\boldsymbol r},{\boldsymbol r}_s^{(p)})$, where $G({\boldsymbol r},{\boldsymbol r}^\prime)$ is Green’s function, which satisfies Eq. (2) with the source term $f(x) = \delta ({\boldsymbol r} - {\boldsymbol r}^\prime)$. Let us subtract the equation for ${u_0}$ from the equation for $u$ [3]:

$$\left\{{\begin{array}{*{20}{l}}- {D_0}\Delta \left({u({\boldsymbol r}) - {u_0}({\boldsymbol r})} \right) + {{\bar\mu}_a}\left({u({\boldsymbol r}) - {u_0}({\boldsymbol r})} \right) = - \delta {\mu _a}({\boldsymbol r})u({\boldsymbol r}),& {\boldsymbol r} \in \Omega ,\\{D_0}{\boldsymbol \nu} \cdot \nabla \left({u({\boldsymbol r}) - {u_0}({\boldsymbol r})} \right) + \frac{1}{\zeta}\left({u({\boldsymbol r}) - {u_0}({\boldsymbol r})} \right) = 0,& {\boldsymbol r} \in \partial \Omega .\end{array}} \right.\\$$
Since the above diffusion equation for $u - {u_0}$ has the source term ${-}\delta {\mu _a}u$, the solution $u({\boldsymbol r})$ satisfies the following identity:
$$u({\boldsymbol r}) = {u_0}({\boldsymbol r}) - \int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime)u({\boldsymbol r}^\prime) {\rm d}{\boldsymbol r}^\prime .$$

The $n$th Born approximation ${u_n}$ is given by ${u_n} = \sum\nolimits_{k = 0}^n {v_k}$, where

$${v_{k + 1}}({\boldsymbol r}) = - \int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime){v_k}({\boldsymbol r}^\prime) {\rm d}{\boldsymbol r}^\prime ,\quad k = 0,1, \ldots ,$$
with ${v_0}({\boldsymbol r}) = {u_0}({\boldsymbol r})$. When the perturbation $\delta {\mu _a}$ is small, the Born series converges and $u = \mathop {\lim}\nolimits_{n \to \infty} {u_n}$. In particular, for non-negative $\delta {\mu _a}$, a concise proof of the convergence is known under the condition $0 \le \delta {\mu _a}({\boldsymbol r}) \le {\bar\mu_a}$ for any ${\boldsymbol r} \in \Omega$ [16]. Let us define the first and second Rytov approximations ${u_R},{u_{R2}}$ as
$${u_R} = {u_0}{e^{{v_1}/{u_0}}},\quad {u_{R2}} = {u_0}{e^{{v_1}/{u_0}}}\exp \left[{\frac{{{v_2}}}{{{u_0}}} - \frac{1}{2}{{\left({\frac{{{v_1}}}{{{u_0}}}} \right)}^2}} \right].$$

When the first Born approximation is compared with the first Rytov approximation, the superiority of the latter has been discussed [3,17,18].

The difference $u - {u_R}$ can be written as

$$\begin{split}{u({\boldsymbol r}) - {u_R}({\boldsymbol r})}&{= {u_0}({\boldsymbol r})\left({1 - {e^{{v_1}({\boldsymbol r})/{u_0}({\boldsymbol r})}}} \right)}\\ &\quad-{ \int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime){u_R}({\boldsymbol r}^\prime) d{\boldsymbol r}^\prime}\\ &\quad- {\int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime)\left({u({\boldsymbol r}^\prime) - {u_R}({\boldsymbol r}^\prime)} \right) d{\boldsymbol r}^\prime .}\end{split}$$

When $\parallel {\mu _a}{\parallel _{{L^\infty}(\Omega)}}$ is small, we have

$$1 - {e^{{v_1}({\boldsymbol r})/{u_0}({\boldsymbol r})}} \approx - \frac{{{v_1}({\boldsymbol r})}}{{{u_0}({\boldsymbol r})}} = \frac{1}{{{u_0}({\boldsymbol r})}}\int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime){u_0}({\boldsymbol r}^\prime) {\rm d}{\boldsymbol r}^\prime ,$$
and
$$\begin{split}&\int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime){u_R}({\boldsymbol r}^\prime) {\rm d}{\boldsymbol r}^\prime \\ &\quad= \int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime){u_0}({\boldsymbol r}^\prime){e^{{v_1}({\boldsymbol r}^\prime)/{u_0}({\boldsymbol r}^\prime)}} {\rm d}{\boldsymbol r}^\prime \\ &\quad\approx \int_\Omega G({\boldsymbol r},{\boldsymbol r}^\prime)\delta {\mu _a}({\boldsymbol r}^\prime){u_0}({\boldsymbol r}^\prime) {\rm d}{\boldsymbol r}^\prime .\end{split}$$

Therefore, for sufficiently small $\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}}$, we obtain

$$\begin{split}\parallel u - {u_R}{\parallel _{{L^\infty}(\Omega)}} &\le {C_{1^\prime}}\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}}\parallel {u_0}{\parallel _{{L^\infty}(\Omega)}}\\ &+ {C_{1^\prime}}\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}}\parallel u - {u_R}{\parallel _{{L^\infty}(\Omega)}},\end{split}$$
where ${C^\prime _1}$ is a positive constant. By moving the last term on the right-hand side to the left-hand side, we arrive at the following inequality for sufficiently small $\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}}$:
$$\parallel u - {u_R}{\parallel _{{L^\infty}(\Omega)}} \le {C_1}\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}}\parallel {u_0}{\parallel _{{L^\infty}(\Omega)}},$$
where ${C_1}$ is a positive constant.

The outgoing light is detected at ${\boldsymbol r}_d^{(p)} \in \partial \Omega$. Let us introduce the data ${\phi ^{(p)}}$ as

$${\phi ^{(p)}} = \ln \frac{{{u_0}({\boldsymbol r}_d^{(p)})}}{{u({\boldsymbol r}_d^{(p)})}}.$$

In the first Rytov approximation, we have ${\phi ^{(p)}} \approx \phi _{R2}^{(p)} \approx \phi _R^{(p)}$, where

$$\begin{split}\phi _R^{(p)} & = - \frac{{{v_1}({\boldsymbol r}_d^{(p)})}}{{{u_0}({\boldsymbol r}_d^{(p)})}},\\ \phi _{R2}^{(p)} & = - \frac{{{v_1}({\boldsymbol r}_d^{(p)})}}{{{u_0}({\boldsymbol r}_d^{(p)})}} + \frac{1}{2}{{\left({\frac{{{v_1}({\boldsymbol r}_d^{(p)})}}{{{u_0}({\boldsymbol r}_d^{(p)})}}} \right)}^2} - \frac{{{v_2}({\boldsymbol r}_d^{(p)})}}{{{u_0}({\boldsymbol r}_d^{(p)})}}.\end{split}$$

Note that ${u_0},{u_R}$ are positive. Using Eq. (14), we have

$$\begin{split}\left| {{e^{- {\phi ^{(p)}}}}} \right| & = \left| {\frac{{u - {u_R}}}{{{u_0}}} + \frac{{{u_R}}}{{{u_0}}}} \right|\\ &\le \frac{{{C_1}\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}}\parallel {u_0}{\parallel _{{L^\infty}(\Omega)}}}}{{{u_0}}} + {e^{{v_1}/{u_0}}}.\end{split}$$

The above inequality implies that ${\phi ^{(p)}} \to 0$ as $\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}} \to 0$.

4. DISCRETIZATION

Let $M$ be an even integer. We introduce a spin variable of discrete values as

$$S({\boldsymbol r}) = 0, \pm 1, \pm 2, \ldots , \pm \frac{M}{2},\quad {\boldsymbol r} \in {\Omega _{{\rm ROI}}}.$$

To explain how $S({\boldsymbol r})$ can be related to $\delta {\mu _a}({\boldsymbol r})$, we will consider the following two cases. In both cases, the two-dimensional space ($d = 2$) is assumed.

A. Single-Spin Model

Let us assume

$$\delta {\mu _a}({\boldsymbol r}) = \eta {f_{\bar a}}(x)\delta (y - {y_0}),$$
where $\eta ,{y_0}$ are given positive constants. In this case, ${\Omega _{{\rm ROI}}}$ is a point $(x,y) = (0,{y_0})$. Here, we assume that ${f_{\bar a}}(x)$ is given by
$${f_{\bar a}}(x) = \left[{{{\bar a}^3} + 3\left({1 + \frac{{\tanh {x^2}}}{{10}}} \right){{\bar a}^2}} \right]\left({1 - \tanh {x^2}} \right),$$
where $\bar a$ is a constant. Thus, $\delta {\mu _a}$ is determined by $\bar a$. This $\bar a$ is the unknown parameter to be reconstructed in the forward data.

When $\bar a$ is sought, we assume that the minimum ${a^{{\rm (min})}}$ and maximum ${a^{{\rm (max})}}$ of candidates of $\bar a$ are known a priori. Let $a$ be a candidate of $\bar a$. We give $a$ as

$$a = {a^{{\rm (min})}} + \frac{{{a^{{\rm (max})}} - {a^{{\rm (min})}}}}{{M + 1}}\left({S + \frac{M}{2} + 1} \right) \in ({a^{{\rm (min})}},{a^{{\rm (max})}}].$$

B. Multi-Spin Model

Next, we suppose that the support of $\delta {\mu _a}$ is unknown, but we know $\delta {\mu _a}$ is zero outside a domain ${\Omega _{{\rm ROI}}}$ ($\mathop {{\rm supp}\,}\delta {\mu _a} \subset {\Omega _{{\rm ROI}}}$).

We divide the region of interest ${\Omega _{{\rm ROI}}}$ into cells ${\omega _i} \subset \Omega$ ($i = 1, \ldots ,N$) such that ${\Omega _{{\rm ROI}}} = \bigcup\nolimits_{i = 1}^N {\omega _i}$ and ${\omega _{{i_1}}} \cap {\omega _{{i_2}}} = \emptyset$ if ${i_1} \ne {i_2}$. Let $|\omega |$ denote the area (volume) of subdomains ${\omega _i}$. We suppose $N \gt {M_{{\rm SD}}}$, which means that the inverse problem is underdetermined. Let us assume $\delta {\mu _a} \in [0,\delta\mu_a^{{\rm (max})}]$ with a constant $\delta\mu_a^{{\rm (max})} \gt 0$. We introduce

$$\delta {\mu _a}({{\boldsymbol r}_i}) = \delta\mu_a^{{\rm (max})}\left({\frac{{{S_i}}}{M} + \frac{1}{2}} \right),\quad {S_i} = 0, \pm 1, \pm 2, \ldots , \pm \frac{M}{2}$$
for $i = 1, \ldots ,N$. Here, ${{\boldsymbol r}_i} \in {\omega _i}$ is a representative point in ${\omega _i}$. Thus, $\delta {\mu _a}$ is discretized in $[0,\delta\mu_a^{{\rm (max})}]$ by $M + 1$ values at each point ${{\boldsymbol r}_i}$.

5. SPIN HAMILTONIAN: CASE 1

We consider the single-spin model of Eq. (19) for the time-dependent diffusion of Eq. (1).

We note that Green’s function for Eq. (1) is given by

$$\begin{split}G({\boldsymbol r},t;{\boldsymbol r}^\prime ,s) & = \frac{{{e^{- {{\bar\mu}_a}c(t - s)}}}}{{4\pi {D_0}(t - s)}}{e^{- \frac{{{{(x - x^\prime)}^2}}}{{4{D_0}c(t - s)}}}}\left[{e^{- \frac{{{{(y - y^\prime)}^2}}}{{4{D_0}c(t - s)}}}} + {e^{- \frac{{{{(y + y^\prime)}^2}}}{{4{D_0}c(t - s)}}}}\right.\\ & \quad- \frac{{\sqrt {4\pi {D_0}c(t - s)}}}{\ell}{e^{- \frac{{{{(y + y^\prime)}^2}}}{{4{D_0}c(t - s)}}}}{e^{{{\left({\frac{{y + y^\prime + 2{D_0}c(t - s)/\ell}}{{\sqrt {4{D_0}c(t - s)}}}} \right)}^2}}}\\&\quad\times\left.{\rm erfc}\left({\frac{{y + y^\prime + 2{D_0}c(t - s)/\ell}}{{\sqrt {4{D_0}c(t - s)}}}} \right)\right]\end{split}$$
for $t \gt s$, and otherwise $G({\boldsymbol r},t;{\boldsymbol r}^\prime ,s) = 0$. Here, $\ell = {D_0}\zeta$ and the complementary error function is given by ${{\rm erfc}} (x) = (2/\sqrt \pi)\int_x^\infty \exp (- {t^2}) {\rm d}t$. Using Green’s function, we obtain
$$\begin{split}{u_0}({\boldsymbol r},t) &= \frac{{{g_0}{e^{- {{\bar{\mu}}_a}ct}}}}{{2\pi {D_0}t}}{e^{- \frac{{{{(x - x_s^i)}^2} + {y^2}}}{{4{D_0}ct}}}}\\&\quad\times\left[{1 - \frac{{\sqrt {\pi {D_0}ct}}}{\ell}{e^{{{\left({\frac{{y + 2{D_0}ct/\ell}}{{\sqrt {4{D_0}ct}}}} \right)}^2}}}{{\rm erfc}}\left({\frac{{y + 2{D_0}ct/\ell}}{{\sqrt {4{D_0}ct}}}} \right)} \right].\end{split}$$

Let us define

$$\begin{split}{g(y,t;y^\prime ,s)}&= {\frac{1}{{4\pi {D_0}(t - s)}}{e^{- \frac{{{{(y + y^\prime)}^2}}}{{4{D_0}c(t - s)}}}}\left[1 + {e^{\frac{{{{(y + y^\prime)}^2} - {{(y - y^\prime)}^2}}}{{4{D_0}c(t - s)}}}}\right.} \\&\quad- {\frac{{\sqrt {4\pi {D_0}c(t - s)}}}{\ell}}{{e^{{{\left({\frac{{y + y^\prime}}{{2\sqrt {{D_0}c(t - s)}}} + \frac{{\sqrt {{D_0}c(t - s)}}}{\ell}} \right)}^2}}}}\\&\quad\times{\left.{\rm erfc}\left({\frac{{y + y^\prime}}{{2\sqrt {{D_0}c(t - s)}}} + \frac{{\sqrt {{D_0}c(t - s)}}}{\ell}} \right)\right].}\end{split}$$

Then, we have

$$\begin{split}&\int_0^t G({\boldsymbol r},t;{\boldsymbol r}^\prime ,s){u_0}({\boldsymbol r}^\prime ,s) {\rm d}s = {e^{- {{\bar\mu}_a}ct}}\int_0^t {e^{- \frac{{{{(x - x^\prime)}^2}}}{{4{D_0}c(t - s)}}}}{e^{- \frac{{{{(x^\prime - x_s^{(p)})}^2}}}{{4{D_0}cs}}}}g(y,t;y^\prime ,s)g(y^\prime ,s;0,0) {\rm d}s.\end{split}$$

Therefore, we obtain

$$\begin{split}{{u_R}({\boldsymbol r},t;{\boldsymbol r}_s^{(p)})}= {{u_0}({\boldsymbol r},t;{\boldsymbol r}_s^{(p)})\exp \left[- \frac{{{e^{- {{\bar\mu}_a}ct}}}}{{{u_0}({\boldsymbol r},t;{\boldsymbol r}_s^i)}}\int_0^\infty \int_{- \infty}^\infty \delta {\mu _a}({\boldsymbol r}^\prime)\left({\int_0^t {e^{- \frac{{{{(x - x^\prime)}^2}}}{{4{D_0}c(t - s)}}}}{e^{- \frac{{{{(x^\prime - x_s^{(p)})}^2}}}{{4{D_0}cs}}}}g(y,t;y^\prime ,s)g(y^\prime ,s;0,0) ds} \right) dx^\prime dy^\prime\right].}\end{split}$$

We will use Eq. (27) to obtain the forward data in our numerical experiment.

By substituting the form of Eq. (19) for $\delta {\mu _a}$ in Eq. (27), we obtain

$$\begin{split}{u_R}({\boldsymbol r}_d^{(p)},t;{\boldsymbol r}_s^{(p)};a) & = u({\boldsymbol r}_d^{(p)},t)\\ & = {u_0}({\boldsymbol r}_d^{(p)},t;{\boldsymbol r}_s^{(p)})\exp \left[- \frac{{\eta {e^{- {{\bar\mu}_a}ct}}}}{{{u_0}({\boldsymbol r}_d^{(p)},t;{\boldsymbol r}_s^{(p)})}}\int_0^t g(0,t;{y_0},s)\right. \left.g({y_0},s;0,0)\left({\int_{- \infty}^\infty {f_a}(x^\prime){e^{- \frac{{{{(x_d^{(p)} - x^\prime)}^2}}}{{4{D_0}c(t - s)}}}}{e^{- \frac{{{{(x^\prime - x_s^{(p)})}^2}}}{{4{D_0}cs}}}} {\rm d}x^\prime} \right) {\rm d}s\right],\end{split}$$
where
$$\begin{split}{u_0} = \frac{{{g_0}{e^{- {{\bar\mu}_a}ct}}}}{{2\pi {D_0}t}}{e^{- \frac{{{{(x_d^{(p)} - x_s^{(p)})}^2}}}{{4{D_0}ct}}}}\left[{1 - \frac{{\sqrt {\pi {D_0}ct}}}{\ell}{e^{{{\left({\frac{{\sqrt {{D_0}ct}}}{\ell}} \right)}^2}}}{\rm erfc}\left({\frac{{\sqrt {{D_0}ct}}}{\ell}} \right)} \right].\end{split}$$

We obtain

$$\begin{split}{\phi _R^{(p)}}&= {\frac{{\eta {e^{- {{\bar\mu}_a}ct}}}}{{{u_0}({\boldsymbol r}_d^{(p)},t;{\boldsymbol r}_s^{(p)})}}\int_0^t g(0,t;{y_0},s)g({y_0},s;0,0)}{\left({\int_{- \infty}^\infty {f_a}(x^\prime){e^{- \frac{{{{(x_d^{(p)} - x^\prime)}^2}}}{{4{D_0}c(t - s)}}}}{e^{- \frac{{{{(x^\prime - x_s^{(p)})}^2}}}{{4{D_0}cs}}}} {\rm d}x^\prime} \right) {\rm d}s.}\end{split}$$

Let us discretize time as $t \to {t_j}$ ($j = 1, \ldots ,{M_t}$). With this discretization, we write $\tilde \phi _R^{(p)}$ instead of $\phi _R^{(p)}$.

In order to use SA instead of the naive Metropolis–Hastings Markov-chain Monte Carlo (MCMC) [10], we treat the cost function as a spin Hamiltonian. We let ${\Phi ^{(p)}}$ denote the experimentally obtained data corresponding to ${\phi ^{(p)}}$. We can solve the inverse problem for our DOT by minimizing the cost function, which is given by

$${\cal H}(S) = \frac{1}{2}\sum\limits_{p = 1}^{{M_{{\rm SD}}}} \sum\limits_{j = 1}^{{M_t}} {\left| {{\Phi ^{(p)}} - \tilde \phi _R^{(p)}} \right|^2} + \alpha |S - {S^{(0)}}|,$$
where $\alpha$ is the regularization parameter, and ${S^{(0)}}$ is an initial guess. We put $\alpha = 0$.

The cost function in Eq. (31) has one local minimum and one global minimum [10]. To see this structure of ${\cal H}$, we introduce

$$\begin{split}h(t) = \frac{1}{t}\exp \left({- \frac{{y_0^2}}{{4Dct}}} \right)\left[{1 - \frac{{\sqrt {\pi Dct}}}{\ell}{e^{{{\left({\frac{{{y_0}}}{{2\sqrt {Dct}}} + \frac{{\sqrt {Dct}}}{\ell}} \right)}^2}}}{\rm erfc}\left({\frac{{{y_0}}}{{2\sqrt {Dct}}} + \frac{{\sqrt {Dct}}}{\ell}} \right)} \right].\end{split}$$

The following form is obtained using Eq. (28). By neglecting noise, we have

$$\begin{split}&u\left({{\boldsymbol r}_d^{(p)},{t_j};{\boldsymbol r}_s^{(p)};\bar a} \right) - u\left({{\boldsymbol r}_d^{(p)},{t_j};{\boldsymbol r}_s^{(p)};a} \right)= \frac{{\eta {e^{- {\mu _{a0}}c{t_j}}}}}{{{{(2\pi D)}^2}}}\int_0^{{t_j}} h({t_j} - s)h(s) \left[{\int_{- \infty}^\infty {d_{\bar a}}(a,x^\prime)\left({1 - \tanh {{x^\prime}^2}} \right){e^{- \frac{{{{(x_d^{(p)} - x^\prime)}^2}}}{{4Dc({t_j} - s)}}}}{e^{- \frac{{{{(x^\prime - x_s^{(p)})}^2}}}{{4Dcs}}}} {\rm d}x^\prime} \right] {\rm d}s,\end{split}$$
where ${d_{\bar a}}(a,x^\prime) = \xi (\bar a,x^\prime) - \xi (a,x^\prime)$ with
$$\xi (a,x^\prime) = {a^2}\left[{a + 3\left({1 + \frac{{\tanh {{x^\prime}^2}}}{{10}}} \right)} \right].$$

For a given $x^\prime $, the function $|{d_{\bar a}}(a;x^\prime {)|^2}$ has the global minimum at $a = \bar a \gt 0$, a local minimum at $a = - 2({1 + \frac{{\tanh {{x^\prime}^2}}}{{10}}})\lt 0$, and a local maximum at $a = 0$. This structure of $|{d_{\bar a}}(a;x^\prime {)|^2}$ implies that iterative methods, which always look for a position that lowers the value of the cost function, do not work when the initial guess ${a_0}$ is negative.

6. SPIN HAMILTONIAN: CASE 2

Here, we consider the multi-spin model of Eq. (22) for the time-independent diffusion of Eq. (2).

Let us define

$${K_{p,i}} = \delta\mu_a^{{\rm (max})}|\omega |\frac{{G({\boldsymbol r}_d^{(p)},{{\boldsymbol r}_i})G({{\boldsymbol r}_i},{\boldsymbol r}_s^{(p)})}}{{G({\boldsymbol r}_d^{(p)},{\boldsymbol r}_s^{(p)})}},$$
where $G$ is Green’s function for Eq. (2) with ${\mu _a} = {\bar\mu_a}$. Let $\tilde \phi _R^{(p)}$, $\tilde \phi _{R2}^{(p)}$ be discretized versions of $\phi _R^{(p)}$, $\phi _{R2}^{(p)}$, respectively. Assuming that ${\omega _i}$ is small, we have
$$\phi _R^{(p)} \approx \tilde \phi _R^{(p)},\quad \phi _{R2}^{(p)} \approx \tilde \phi _{R2}^{(p)},$$
where
$$\tilde \phi _R^{(p)} = \sum\limits_{i = 1}^N {K_{p,i}}\left({\frac{{{S_i}}}{M} + \frac{1}{2}} \right),$$
and
$$\begin{split}{\tilde \phi _{R2}^{(p)}}&{= \sum\limits_{i = 1}^N {K_{p,i}}\left({\frac{{{S_i}}}{M} + \frac{1}{2}} \right)}\\ &+{ \frac{1}{2}\sum\limits_{{i_1} = 1}^N \sum\limits_{{i_2} = 1}^N {K_{p,{i_1}}}{K_{p,{i_2}}}\left({\frac{{{S_{{i_1}}}}}{M} + \frac{1}{2}} \right)\left({\frac{{{S_{{i_2}}}}}{M} + \frac{1}{2}} \right)}\\ &\quad-{ |\omega {|^2}\sum\limits_{{i_1} = 1}^N \sum\limits_{{i_2} = 1}^N {{\left[{G({\boldsymbol r}_d^{(p)},{\boldsymbol r}_s^{(p)})} \right]}^{- 1}}}\\ &\quad\times{ G({\boldsymbol r}_d^{(p)},{{\boldsymbol r}_{{i_1}}})\delta {\mu _a}({{\boldsymbol r}_{{i_1}}})G({{\boldsymbol r}_{{i_1}}},{{\boldsymbol r}_{{i_2}}})\delta {\mu _a}({{\boldsymbol r}_{{i_2}}})G({{\boldsymbol r}_{{i_2}}},{\boldsymbol r}_s^{(p)}).}\end{split}$$

Let ${\boldsymbol S} = ({S_i})$ ($i = 1, \ldots ,N$) denote the spin configuration. We will look for ${{\boldsymbol S}^*}$, which minimizes the following cost function $\Psi ({\boldsymbol S})$:

$$\Psi ({\boldsymbol S}) = \frac{1}{2}\sum\limits_{p = 1}^{{M_{{\rm SD}}}} {\left| {{\Phi ^{(p)}} - \tilde \phi _{R2}^{(p)}} \right|^2} + \alpha \sum\limits_{i = 1}^N |{S_i} - S_i^{(0)}|,$$
where $\alpha \gt 0$ is the regularization parameter, and ${{\textbf{S}}^{(0)}}=(S_{i}^{(0)})\in {{\mathbb{R}}^{N}}$ is an initial guess, which we set as $S_i^{(0)} = - M/2$ ($i = 1, \ldots ,N$).

To compare terms of order ${(\delta\mu_a^{{\rm (max)}})^2}$, let us introduce

$$\begin{split}{{g_1}}&= {{K_{p,{i_1}}}{K_{p,{i_2}}},}\\{{g_2}}&={ {g_1} - 2{{\left({\delta\mu_a^{{\rm (max})}} \right)}^2}|\omega {|^2}\frac{{G({\boldsymbol r}_d^{(p)},{{\boldsymbol r}_{{i_1}}})G({{\boldsymbol r}_{{i_1}}},{{\boldsymbol r}_{{i_2}}})G({{\boldsymbol r}_{{i_2}}},{\boldsymbol r}_s^{(p)})}}{{G({\boldsymbol r}_d^{(p)},{\boldsymbol r}_s^{(p)})}}.}\end{split}$$

We have

$$\begin{split}{\left| {{\phi ^{(p)}}{g_2}} \right|}&{= |{g_1}|\left| {{\phi ^{(p)}}\frac{{{g_2}}}{{{g_1}}}} \right|}\\ &{= |{g_1}|\left| {{\phi ^{(p)}}} \right|\left| {1 - 2\frac{{G({\boldsymbol r}_d^{(p)},{\boldsymbol r}_s^{(p)})G({{\boldsymbol r}_{{i_1}}},{{\boldsymbol r}_{{i_2}}})}}{{G({\boldsymbol r}_d^{(p)},{{\boldsymbol r}_{{i_2}}})G({{\boldsymbol r}_{{i_1}}},{\boldsymbol r}_s^{(p)})}}} \right|.}\end{split}$$

Recall that $|{\phi ^{(p)}}| \to 0$ as $\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}} \to 0$. Thus, for arbitrary $p,{i_1},{i_2}$ ($p = 1, \ldots ,{M_{{\rm SD}}}$, ${i_1},{i_2} = 1, \ldots ,N$),

$$|{g_1}| \ge |{\phi ^{(p)}}{g_2}|$$
when $\parallel \delta {\mu _a}{\parallel _{{L^\infty}(\Omega)}}$ is sufficiently small.

By neglecting ${g_2}$, we have $\Psi = {\cal H}({\boldsymbol S}) + {\rm const}.$, where only the first term ${\cal H}({\boldsymbol S})$ on the right-hand side depends on ${S_i}$. The Hamiltonian ${\cal H}({\boldsymbol S})$ is given by

$${\cal H}({\boldsymbol S}) = - \sum\limits_{i = 1}^N \sum\limits_{j = 1}^N {J_{\textit{ij}}}{S_i}{S_j} - \sum\limits_{i = 1}^N {h_i}{S_i},$$
where
$$\begin{split}{{J_{\textit{ij}}}}&= {\frac{{- 1}}{{2{M^2}}}\sum\limits_{p = 1}^{{M_{{\rm SD}}}} {K_{p,i}}{K_{p,j}},}\\{{h_i}}&= {M\sum\limits_{j = 1}^N {J_{\textit{ij}}} + \frac{1}{M}\left({\sum\limits_{p = 1}^{{M_{{\rm SD}}}} {\Phi ^{(p)}}{K_{p,i}} - \alpha} \right).}\end{split}$$
The spin interactions are symmetric: ${J_{\textit{ij}}} = {J_{\textit{ji}}}$. We note that the modulus $|{h_i}|$ of the magnetic field becomes large when the hyperparameter $\alpha$ in the regularization term is large. In this case, the spin Hamiltonian has a unique ground state (${{\boldsymbol S}^*} = {{\boldsymbol S}^{(0)}}$) for sufficiently large $\alpha$.

Thus, our optical tomography is reformulated as the problem of searching the ground state of ${\cal H}({\boldsymbol S})$, i.e., the spin configuration ${{\boldsymbol S}^*} = {{\rm argmin}}_S {\cal H}({\boldsymbol S})$.

7. SIMULATED ANNEALING

The algorithm of SA is described below. Although we assume Case 2 of multiple spins, the single spin in Case 1 can be similarly implemented.

The Hamiltonian ${\cal H}({\boldsymbol S})$ in Eq. (43) can be seen as a Hamiltonian of spins that interact with each other in a solid state material. In a real physical spin system, spins fluctuate depending on temperature, and the spin system reaches its lowest energy state if temperature gradually decreases. In SA, we seek the optimal solution by mimicking this physical process with thermal fluctuations. To this end, we introduce the parameter $T$, which corresponds to temperature. We will find the ground state ${{\boldsymbol S}^*}$ of ${\cal H}({\boldsymbol S})$ in Eq. (43) by gradually decreasing $T$ from ${T_{{\rm high}}}$ to ${T_{{\rm low}}}$.

The partition function $Z$ is given by

$$Z = \sum\limits_{\{{S_i}\}} {e^{- \beta {\cal H}({\boldsymbol S})}},$$
where $\beta = 1/T$ is the inverse temperature. Here, we used the notation $\sum\nolimits_{\{{S_i}\}} = \sum\nolimits_{{S_1}} \cdots \sum\nolimits_{{S_N}}$. We treat $\{{S_i}\} = \{{S_1},{S_2}, \ldots ,{S_N}\}$ as random variables. We give the probability density function $\pi ({\boldsymbol S})$ as the Boltzmann distribution given by
$$\pi ({\boldsymbol S}) = \frac{{{e^{- \beta {\cal H}({\boldsymbol S})}}}}{Z}.$$

We see that $\pi ({\boldsymbol S})$ is large if ${\cal H}({\boldsymbol S})$ is small. When the temperature is sufficiently low, i.e., $\beta$ is large, $\pi ({{\boldsymbol S}^*})$ becomes significantly larger than the probability of other configurations.

Although calculating the denominator $Z$ on the right-hand side of Eq. (46) is difficult, we can find ${{\boldsymbol S}^*}$ by using the Metropolis algorithm [15,19]. The proposal distribution $q({S^\prime _i}|{\boldsymbol S})$ is given for each spin, say the $i$th spin, and the value of ${S^\prime _i}$ is generated with equal probability. We note that

$$\begin{split}\frac{{\pi ({\boldsymbol S}^\prime)}}{{\pi ({\boldsymbol S})}} & = {e^{\beta ({\cal H}({\boldsymbol S}) - {\cal H}({\boldsymbol S}^\prime))}}\\ & = \exp \left[{\beta \left({{J_{\textit{ii}}}\left({{{S^\prime_i}}^2 - S_i^2} \right) + \left({2\sum\limits_{\begin{array}{*{20}{c}}{j = 1}\\{j \ne i}\end{array}}^N {J_{\textit{ij}}}{S_j} + {h_i}} \right)\left({{{S^\prime_i}} - {S_i}} \right)} \right)} \right].\end{split}$$

While the naive Metropolis–Hastings algorithm does not work in high dimensions [20,21], the Metropolis algorithm or single-component Metropolis–Hastings algorithm can handle many unknown parameters [22]. Although we use the standard Metropolis algorithm as the first attempt of the single-component Metropolis–Hastings algorithm for DOT, different algorithms with faster relaxation have been developed for spin systems [15]. Our algorithm is summarized as the following six steps:

  • 1. Start with a small $\beta = 1/{T_{{\rm high}}} \gt 0$. Give $({S_i})$ randomly as an initial guess. Then, set $i = 1$.
  • 2. Compute ${h_{{\rm eff},i}} = 2\sum\nolimits_{j = 1\;(j \ne i)}^N {J_{\textit{ij}}}{S_j} + {h_i}$.
  • 3. Calculate $w = - \beta [{h_{{\rm eff},i}}({S^\prime _i} - {S_i}) + {J_{\textit{ii}}}({S^\prime _i}^2 - S_i^2)]$, where ${S^\prime _i} \sim q(\cdot |{\boldsymbol S})$ is randomly chosen.
  • 4. Set ${S_i} = {S^\prime _i}$ if $w \le 0$. Otherwise, put ${S_i} = {S^\prime _i}$ with probability ${e^{- w}}$.
  • 5. Set $i = 1$ if $i = N$. Otherwise, set $i = i + 1$. Return to Step 2. After repeating several loops from Step 2 to Step 5 until the initial large fluctuation ceases, proceed to Step 6.
  • 6. Decrease temperature and go to Step 2. If the temperature reaches ${T_{{\rm low}}}$, finish the iteration.

In this paper, we decrease $T$ as

$$T - {10^{{\rm int}(\mathop {\log}\nolimits_{10} T) - 2}} \to T.$$

8. NUMERICAL TEST

A. Single Spin

We take the unit of length and unit of time to be millimeters (mm) and picoseconds (ps), respectively. On the $x$ axis, we place two sources at $(x,y) = (- 20,0)$, $(20,0)$ and three detectors at $(x,y) = (- 40,0)$, (0,0), (40,0). As a result, we have ${M_{{\rm SD}}} = 6$. We set ${D_0} = 0.33$, ${\mu _a} = 0.02$, and $\mathfrak{n}=1.37$. Suppose that there is absorption inhomogeneity at depth five. For $\delta {\mu _a}$, we put $\eta = 300/c$, ${y_0} = 5$, and

$$\bar a = 1.5.$$

When we prepare ${\Phi ^{(p)}}$, we added $3\%$ Gaussian noise. We set

$${t_j} = j{\Delta _t}\quad (j = 1, \ldots ,{M_t}),\quad {\Delta _t} = 5,\quad {M_t} = 500.$$

Furthermore, we set $M = 512$, ${a^{{\rm (min})}} = - 3$, and ${a^{{\rm (max})}} = 3$.

It is known that iterative methods cannot reach $\bar a = 1.5$ when the initial value is negative [10]. We set the initial value ${a_0}$ of $a$ as

$${a_0} = - 0.01.$$

In Fig. 1, we compare the SA developed in this paper with the Levenberg–Marquardt algorithm [2325], which is one of iterative methods. The Levenberg–Marquardt algorithm fails to arrive at the correct answer, whereas SA converges to the correct value. In our simulation, converged values for the Levenberg–Marquardt algorithm and SA are ${-}2.05$ and 1.68, respectively. As is mentioned in the end of Section 5, the cost function has a local minimum at a negative value of $a$. The calculation of the Levenberg–Marquardt algorithm is trapped by this local minimum. On the other hand, $a$ goes back and forth between the local minimum and global minimum in SA and then falls in the global minimum as temperature decreases.

 figure: Fig. 1.

Fig. 1. (top) The value $a$ is plotted against Monte Carlo steps for the simulated annealing in Section 7. (bottom) The reconstruction is done with the Levenberg–Marquardt algorithm.

Download Full Size | PPT Slide | PDF

 figure: Fig. 2.

Fig. 2. (left) One absorber at the depth 10 mm. (right) One absorber at a deeper position at $y = 15\;{\rm mm} $.

Download Full Size | PPT Slide | PDF

B. Multiple Spins

Let us consider DOT in the half space $\Omega$, i.e., $\Omega = \{ \textbf{r}\in {{\mathbb{R}}^{2}}; -\infty \lt x \lt \infty ,\;0\lt y\lt \infty \}$. The boundary $\partial \Omega$ is the $x$ axis. Measurements are performed on the $x$ axis: ${\boldsymbol r}_s^{(p)} = (x_s^{(p)},0^{+}), {\boldsymbol r}_d^{(p)} = (x_d^{(p)},0)$.

 figure: Fig. 3.

Fig. 3. Two absorbers with separation of 20 mm reconstructed by the proposed method.

Download Full Size | PPT Slide | PDF

 figure: Fig. 4.

Fig. 4. Two absorbers with separation of 20 mm. (left) The reconstruction by the truncated SVD with 52 singular values. (right) The reconstruction by the truncated SVD with 80 singular values.

Download Full Size | PPT Slide | PDF

In this case, Green’s function is explicitly given by
$$\begin{split}{G({\boldsymbol r},{\boldsymbol r}^\prime)}&{= \frac{1}{{2\pi {D_0}}}\int_0^\infty \frac{{\cos (q(x - x^\prime))}}{\lambda}}\\ &\quad\times {\left({{e^{- \lambda |y - y^\prime |}} + \frac{{\ell \lambda - 1}}{{\ell \lambda + 1}}{e^{- \lambda (y + y^\prime)}}} \right) {\rm d}q,}\end{split}$$
where $\ell = \zeta {D_0}$ and $\lambda = \lambda (q) = \sqrt {\frac{{{{\bar\mu}_a}}}{{{D_0}}} + {q^2}}$. Equation (52) implies $\parallel G({\boldsymbol r}, \cdot){\parallel _{{L^1}(\Omega)}} \lt \infty$ and $\parallel G({\boldsymbol r}, \cdot){\parallel _{{L^\infty}(\Omega)}} \lt \infty$ for any ${\boldsymbol r} \in \Omega$. The integrand in Eq. (52) decays slowly when $y,y^\prime $ are small and is oscillatory. The numerical integration in Eq. (52) can be done by the double-exponential formula [26].

We use ${M_{{\rm SD}}} = 240$ source-detector pairs from 16 sources and 15 detectors:

$$\begin{split}{x_s^{(p)}} &= {\pm 2, \pm 6, \ldots , \pm 30 \;{\rm mm},}\\{x_d^{(p)}} &= {0, \pm 4, \pm 8, \ldots , \pm 28 \;{\rm mm}.}\end{split}$$

In addition to the refractive index $\mathfrak{n}=1.37$, we set

$${\bar\mu_a} = 0.02\;{{\rm mm}^{- 1}},\quad {D_0} = 0.33\;{{\rm mm}^{- 1}}.$$

We assume absorption inhomogeneity of the shape of disks in the medium. Inside the disks, we set

$$\delta {\mu _a} = 0.2\;{{\rm mm}^{- 1}}.$$

We compute the forward data ${\Phi ^{(p)}}$ ($p = 1, \ldots ,{M_{{\rm SD}}}$) by the finite-difference scheme with the Gauss–Seidel method. We added 3% Gaussian noise to $u({\boldsymbol r}_d^{(p)})$ and ${u_0}({\boldsymbol r}_d^{(p)})$. Assuming the diffuse surface reflection, we have [27]

$$\zeta = 2\frac{{1 + {r_d}}}{{1 - {r_d}}},$$
where
$${{r}_{d}}=-1.4399{{\mathfrak{n}}^{-2}}+0.7099{{\mathfrak{n}}^{-1}}+0.6681+0.0636\mathfrak{n}.$$

We take $(2{N_x} + 1){N_y} = 1830$ cells (${N_x} = 30$, ${N_y} = 30$) of the area $|\omega | = {h^2}$ ($h = 1\;{\rm mm} $). Moreover, the following parameter values were used for SA:

$$\alpha = 0.01,\quad {T_{{\rm high}}}{= 10^{- 5}},\quad {T_{{\rm low}}}{= 10^{- 10}},\quad M = 256.$$

Figures 2 and 3 show reconstructed images by the proposed method. In Fig. 2, a disk target of radius 2.5 mm was placed at the depth of 10 mm (left) and 15 mm (right), that is, the $x$ and $y$ coordinates of the disk center are $x = 0\;{\rm mm} $ and $y = 10\;{\rm mm} $ or 15 mm. In general, the reconstruction at deeper positions is more difficult because the signal becomes noisy, and the relative contribution of the regularization term in the cost function becomes larger [12]. As is expected, the reconstruction at a deeper position is more difficult.

Next, we consider two disks of radius of 2.5 mm in Fig. 3. The centers of the disks are placed at $(x,y) = (\pm 10 \;{\rm mm}, 10 \;{\rm mm})$. Figure 3 shows the reconstruction by the proposed scheme using the same parameters described above for Fig. 2. For comparison, we also used the truncated SVD. Reconstructions from the truncated SVD are presented in Fig. 4 with 52 largest singular values (left) and 80 largest singular values (right), respectively. In the reconstructed images by the truncated SVD, the separation between the two targets is less sharp.

In each numerical calculation, there are $5 \times {10^5}$ unknown parameter values. The computation time to find ${{\boldsymbol S}^*}$ was about 3 min on a laptop computer (MacBook Pro with 2.3 GHz Intel Core i5 and 8 GB memory).

9. CONCLUDING REMARKS

We have developed a stochastic approach of DOT using the Metropolis algorithm of the Monte Carlo method. Due to the fact that the inverse problem for DOT is severely ill-posed, the reconstructed values are always blurred at least to some extent. This motivated us to introduce spins, which have discrete values. The effect of the discretization of the parameter is subtle for the single-spin model but becomes important when the number of spins increases in the multi-spin model because the discretization reduces the search space.

For linearized inverse problems, the truncated SVD is often used [3,4]. Although the truncated SVD has a filtering property similar to the Tikhonov regularization, the penalty term of the proposed method is not restricted to the 2-norm. In this paper, the 1-norm was used. This might explain the difference of the quality of reconstruction in Fig. 4. Another advantage of our method is that it is feasible to use a priori information and set the lower and upper bounds of $\delta {\mu _a}$, whereas the reconstructed $\delta {\mu _a}$ by the truncated SVD is unbounded. This feature is reflected in the difference between the proposed method and SVD in Fig. 4. As is seen for the single-spin model, the proposed method can treat both linear and nonlinear inverse problems.

The single-spin model provides a nonlinear inverse problem. Although it is a linearized inverse problem, the multi-spin model demonstrates that a large number of unknowns can be treated by our numerical scheme. It is a natural next step to consider nonlinear inverse problems with multiple spins. For this, nonlinear terms such as ${u_{R2}}$ in the Rytov series have to be taken into account. When the landscape has a more complicated structure, it is important to decrease temperature gradually.

To achieve fast convergence, it is crucial to bring the spin system efficiently to the thermal equilibrium state. Various Monte Carlo methods for spin systems have been developed to reduce the computation time. These techniques may be implemented in our approach, which solves the inverse problem of DOT from the viewpoint of statistical mechanics of spin systems. Existing Monte Carlo methods for spin systems include the following algorithms. The replica exchange method prepares multiple copies of a spin system to let the system at low temperatures escape from local minima [28]. To minimize the average rejection rate, the MCMC algorithm by Suwa and Todo breaks the detailed balance condition while preserving the balance condition [29]. Quantum annealing makes use of quantum fluctuations of spins whereas SA uses thermal fluctuations [30]. Quantum annealing also has potential to provide an efficient algorithm for DOT.

Funding

National Natural Science Foundation of China (11971121); Japan Society for the Promotion of Science KAKENHI (JP17K05572, JP18K03438, JP16K05418); Hamamatsu University School of Medicine.

Acknowledgment

Y. J. is supported by the National Natural Science Foundation of China (No. 11971121). M. M. is supported by JSPS KAKENHI grant numbers JP17K05572 and JP18K03438 and by HUSM Grant-in-Aid. N. T. is supported by JSPS KAKENHI grant number JP16K05418.

Disclosures

The authors declare no conflicts of interest.

Data Availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

REFERENCES

1. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef]  

2. Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev. 21, 185–205 (2014). [CrossRef]  

3. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999). [CrossRef]  

4. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009). [CrossRef]  

5. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E 70, 056616 (2004). [CrossRef]  

6. R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005). [CrossRef]  

7. H. Zhao and R. J. Cooper, “Review of recent progress toward a fiberless, whole-scalp diffuse optical tomography system,” Neurophotonics 5, 011012 (2017). [CrossRef]  

8. Q. Zhu and S. Poplack, “A review of optical breast imaging: multi-modality systems for breast cancer diagnosis,” Eur. J. Radiol. 129, 109067 (2020). [CrossRef]  

9. M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999). [CrossRef]  

10. Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019). [CrossRef]  

11. S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006). [CrossRef]  

12. T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue, and M. Sato, “Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography,” Opt. Express 20, 20427–20446 (2012). [CrossRef]  

13. G. Bal, I. Langmore, and Y. Marzouk, “Bayesian inverse problems with Monte Carlo forward models,” Inverse Probl. Imaging 7, 81–105 (2013). [CrossRef]  

14. I. Langmore, A. B. Davis, and G. Bal, “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. 51, 2903–2919 (2013). [CrossRef]  

15. D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University, 2000).

16. V. A. Markel and J. C. Schotland, “On the convergence of the Born series in optical tomography with diffuse light,” Inverse Probl. 23, 1445–1465 (2007). [CrossRef]  

17. J. B. Keller, “Accuracy and validity of the Born and Rytov approximations,” J. Opt. Soc. Am. 59, 1003–1004 (1969). [CrossRef]  

18. E. Kirkinis, “Renormalization group interpretation of the Born and Rytov approximations,” J. Opt. Soc. Am. A 25, 2499–2508 (2008). [CrossRef]  

19. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953). [CrossRef]  

20. S.-K. Au and J. L. Beck, “Estimation of small failure probabilities in high dimensions by subset simulation,” Prob. Eng. Mech. 16, 263–277 (2001). [CrossRef]  

21. L. S. Katafygiotis and K. M. Zuev, “Geometric insight into the challenges of solving high-dimensional reliability problems,” Prob. Eng. Mech. 23, 208–218 (2008). [CrossRef]  

22. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov Chain Monte Carlo in Practice (Chapman and Hall, 1996).

23. R. Fletcher, “A modified Marquardt subroutine for nonlinear least squares,” Report AERE-R 6799 (The Atomic Energy Research Establishment, Harwell, 1971).

24. K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quart. Appl. Math. 2, 164–168 (1944). [CrossRef]  

25. D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963). [CrossRef]  

26. T. Ooura and M. Mori, “The double exponential formula for oscillatory functions over the half infinite interval,” J. Comput. Appl. Math. 38, 353–360 (1991). [CrossRef]  

27. W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials (Academic, 1979).

28. K. Hukushima and K. Nemoto, “Exchange Monte Carlo method and application to spin glass simulations,” J. Phys. Soc. Jpn. 65, 1604–1608 (1996). [CrossRef]  

29. H. Suwa and S. Todo, “Markov chain Monte Carlo method without detailed balance,” Phys. Rev. Lett. 105, 120603 (2010). [CrossRef]  

30. T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E 58, 5355–5363 (1998). [CrossRef]  

References

  • View by:
  • |
  • |
  • |

  1. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).
    [Crossref]
  2. Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev. 21, 185–205 (2014).
    [Crossref]
  3. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).
    [Crossref]
  4. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009).
    [Crossref]
  5. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E 70, 056616 (2004).
    [Crossref]
  6. R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
    [Crossref]
  7. H. Zhao and R. J. Cooper, “Review of recent progress toward a fiberless, whole-scalp diffuse optical tomography system,” Neurophotonics 5, 011012 (2017).
    [Crossref]
  8. Q. Zhu and S. Poplack, “A review of optical breast imaging: multi-modality systems for breast cancer diagnosis,” Eur. J. Radiol. 129, 109067 (2020).
    [Crossref]
  9. M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999).
    [Crossref]
  10. Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019).
    [Crossref]
  11. S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
    [Crossref]
  12. T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue, and M. Sato, “Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography,” Opt. Express 20, 20427–20446 (2012).
    [Crossref]
  13. G. Bal, I. Langmore, and Y. Marzouk, “Bayesian inverse problems with Monte Carlo forward models,” Inverse Probl. Imaging 7, 81–105 (2013).
    [Crossref]
  14. I. Langmore, A. B. Davis, and G. Bal, “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. 51, 2903–2919 (2013).
    [Crossref]
  15. D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University, 2000).
  16. V. A. Markel and J. C. Schotland, “On the convergence of the Born series in optical tomography with diffuse light,” Inverse Probl. 23, 1445–1465 (2007).
    [Crossref]
  17. J. B. Keller, “Accuracy and validity of the Born and Rytov approximations,” J. Opt. Soc. Am. 59, 1003–1004 (1969).
    [Crossref]
  18. E. Kirkinis, “Renormalization group interpretation of the Born and Rytov approximations,” J. Opt. Soc. Am. A 25, 2499–2508 (2008).
    [Crossref]
  19. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
    [Crossref]
  20. S.-K. Au and J. L. Beck, “Estimation of small failure probabilities in high dimensions by subset simulation,” Prob. Eng. Mech. 16, 263–277 (2001).
    [Crossref]
  21. L. S. Katafygiotis and K. M. Zuev, “Geometric insight into the challenges of solving high-dimensional reliability problems,” Prob. Eng. Mech. 23, 208–218 (2008).
    [Crossref]
  22. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov Chain Monte Carlo in Practice (Chapman and Hall, 1996).
  23. R. Fletcher, “A modified Marquardt subroutine for nonlinear least squares,” Report AERE-R 6799 (The Atomic Energy Research Establishment, Harwell, 1971).
  24. K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quart. Appl. Math. 2, 164–168 (1944).
    [Crossref]
  25. D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963).
    [Crossref]
  26. T. Ooura and M. Mori, “The double exponential formula for oscillatory functions over the half infinite interval,” J. Comput. Appl. Math. 38, 353–360 (1991).
    [Crossref]
  27. W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials (Academic, 1979).
  28. K. Hukushima and K. Nemoto, “Exchange Monte Carlo method and application to spin glass simulations,” J. Phys. Soc. Jpn. 65, 1604–1608 (1996).
    [Crossref]
  29. H. Suwa and S. Todo, “Markov chain Monte Carlo method without detailed balance,” Phys. Rev. Lett. 105, 120603 (2010).
    [Crossref]
  30. T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E 58, 5355–5363 (1998).
    [Crossref]

2020 (1)

Q. Zhu and S. Poplack, “A review of optical breast imaging: multi-modality systems for breast cancer diagnosis,” Eur. J. Radiol. 129, 109067 (2020).
[Crossref]

2019 (1)

Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019).
[Crossref]

2017 (1)

H. Zhao and R. J. Cooper, “Review of recent progress toward a fiberless, whole-scalp diffuse optical tomography system,” Neurophotonics 5, 011012 (2017).
[Crossref]

2014 (1)

Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev. 21, 185–205 (2014).
[Crossref]

2013 (2)

G. Bal, I. Langmore, and Y. Marzouk, “Bayesian inverse problems with Monte Carlo forward models,” Inverse Probl. Imaging 7, 81–105 (2013).
[Crossref]

I. Langmore, A. B. Davis, and G. Bal, “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. 51, 2903–2919 (2013).
[Crossref]

2012 (1)

2010 (1)

H. Suwa and S. Todo, “Markov chain Monte Carlo method without detailed balance,” Phys. Rev. Lett. 105, 120603 (2010).
[Crossref]

2009 (1)

S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009).
[Crossref]

2008 (2)

E. Kirkinis, “Renormalization group interpretation of the Born and Rytov approximations,” J. Opt. Soc. Am. A 25, 2499–2508 (2008).
[Crossref]

L. S. Katafygiotis and K. M. Zuev, “Geometric insight into the challenges of solving high-dimensional reliability problems,” Prob. Eng. Mech. 23, 208–218 (2008).
[Crossref]

2007 (1)

V. A. Markel and J. C. Schotland, “On the convergence of the Born series in optical tomography with diffuse light,” Inverse Probl. 23, 1445–1465 (2007).
[Crossref]

2006 (1)

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

2005 (2)

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).
[Crossref]

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

2004 (1)

V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E 70, 056616 (2004).
[Crossref]

2001 (1)

S.-K. Au and J. L. Beck, “Estimation of small failure probabilities in high dimensions by subset simulation,” Prob. Eng. Mech. 16, 263–277 (2001).
[Crossref]

1999 (2)

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).
[Crossref]

M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999).
[Crossref]

1998 (1)

T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E 58, 5355–5363 (1998).
[Crossref]

1996 (1)

K. Hukushima and K. Nemoto, “Exchange Monte Carlo method and application to spin glass simulations,” J. Phys. Soc. Jpn. 65, 1604–1608 (1996).
[Crossref]

1991 (1)

T. Ooura and M. Mori, “The double exponential formula for oscillatory functions over the half infinite interval,” J. Comput. Appl. Math. 38, 353–360 (1991).
[Crossref]

1969 (1)

1963 (1)

D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963).
[Crossref]

1953 (1)

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
[Crossref]

1944 (1)

K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quart. Appl. Math. 2, 164–168 (1944).
[Crossref]

Amita, T.

Arridge, S. R.

S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009).
[Crossref]

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).
[Crossref]

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).
[Crossref]

M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999).
[Crossref]

Au, S.-K.

S.-K. Au and J. L. Beck, “Estimation of small failure probabilities in high dimensions by subset simulation,” Prob. Eng. Mech. 16, 263–277 (2001).
[Crossref]

Bal, G.

G. Bal, I. Langmore, and Y. Marzouk, “Bayesian inverse problems with Monte Carlo forward models,” Inverse Probl. Imaging 7, 81–105 (2013).
[Crossref]

I. Langmore, A. B. Davis, and G. Bal, “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. 51, 2903–2919 (2013).
[Crossref]

Beck, J. L.

S.-K. Au and J. L. Beck, “Estimation of small failure probabilities in high dimensions by subset simulation,” Prob. Eng. Mech. 16, 263–277 (2001).
[Crossref]

Binder, K.

D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University, 2000).

Chance, B.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Choe, R.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Cooper, R. J.

H. Zhao and R. J. Cooper, “Review of recent progress toward a fiberless, whole-scalp diffuse optical tomography system,” Neurophotonics 5, 011012 (2017).
[Crossref]

Corlu, A.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Czerniecki, B. J.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Davis, A. B.

I. Langmore, A. B. Davis, and G. Bal, “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. 51, 2903–2919 (2013).
[Crossref]

DeMichele, A.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Durduran, T.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Egan, W. G.

W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials (Academic, 1979).

Fletcher, R.

R. Fletcher, “A modified Marquardt subroutine for nonlinear least squares,” Report AERE-R 6799 (The Atomic Energy Research Establishment, Harwell, 1971).

Fraker, D. L.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Gibson, A. P.

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).
[Crossref]

Gilks, W. R.

W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov Chain Monte Carlo in Practice (Chapman and Hall, 1996).

Grosicka-Koptyra, M.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Hebden, J. C.

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).
[Crossref]

Hilgeman, T. W.

W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials (Academic, 1979).

Hiroe, N.

Hoshi, Y.

Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019).
[Crossref]

Hukushima, K.

K. Hukushima and K. Nemoto, “Exchange Monte Carlo method and application to spin glass simulations,” J. Phys. Soc. Jpn. 65, 1604–1608 (1996).
[Crossref]

Inoue, Y.

Jiang, Y.

Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019).
[Crossref]

Kadowaki, T.

T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E 58, 5355–5363 (1998).
[Crossref]

Kaipio, J. P.

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

Katafygiotis, L. S.

L. S. Katafygiotis and K. M. Zuev, “Geometric insight into the challenges of solving high-dimensional reliability problems,” Prob. Eng. Mech. 23, 208–218 (2008).
[Crossref]

Keller, J. B.

Kirkinis, E.

Kolehmainen, V.

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

Konecky, S. D.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Kosaka, T.

Landau, D. P.

D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University, 2000).

Langmore, I.

I. Langmore, A. B. Davis, and G. Bal, “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. 51, 2903–2919 (2013).
[Crossref]

G. Bal, I. Langmore, and Y. Marzouk, “Bayesian inverse problems with Monte Carlo forward models,” Inverse Probl. Imaging 7, 81–105 (2013).
[Crossref]

Lee, K.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Levenberg, K.

K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quart. Appl. Math. 2, 164–168 (1944).
[Crossref]

Machida, M.

Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019).
[Crossref]

Markel, V. A.

V. A. Markel and J. C. Schotland, “On the convergence of the Born series in optical tomography with diffuse light,” Inverse Probl. 23, 1445–1465 (2007).
[Crossref]

V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E 70, 056616 (2004).
[Crossref]

Marquardt, D. W.

D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963).
[Crossref]

Marzouk, Y.

G. Bal, I. Langmore, and Y. Marzouk, “Bayesian inverse problems with Monte Carlo forward models,” Inverse Probl. Imaging 7, 81–105 (2013).
[Crossref]

Metropolis, N.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
[Crossref]

Mori, M.

T. Ooura and M. Mori, “The double exponential formula for oscillatory functions over the half infinite interval,” J. Comput. Appl. Math. 38, 353–360 (1991).
[Crossref]

Nakamura, G.

Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019).
[Crossref]

Nemoto, K.

K. Hukushima and K. Nemoto, “Exchange Monte Carlo method and application to spin glass simulations,” J. Phys. Soc. Jpn. 65, 1604–1608 (1996).
[Crossref]

Nishimori, H.

T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E 58, 5355–5363 (1998).
[Crossref]

Okawa, S.

Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev. 21, 185–205 (2014).
[Crossref]

Ooura, T.

T. Ooura and M. Mori, “The double exponential formula for oscillatory functions over the half infinite interval,” J. Comput. Appl. Math. 38, 353–360 (1991).
[Crossref]

Poplack, S.

Q. Zhu and S. Poplack, “A review of optical breast imaging: multi-modality systems for breast cancer diagnosis,” Eur. J. Radiol. 129, 109067 (2020).
[Crossref]

Richardson, S.

W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov Chain Monte Carlo in Practice (Chapman and Hall, 1996).

Rosen, M. A.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Rosenbluth, A. W.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
[Crossref]

Rosenbluth, M. N.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
[Crossref]

Sato, M.

Schotland, J. C.

S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009).
[Crossref]

V. A. Markel and J. C. Schotland, “On the convergence of the Born series in optical tomography with diffuse light,” Inverse Probl. 23, 1445–1465 (2007).
[Crossref]

V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E 70, 056616 (2004).
[Crossref]

Schweiger, M.

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999).
[Crossref]

Shimokawa, T.

Somersalo, E.

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

Spiegelhalter, D. J.

W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov Chain Monte Carlo in Practice (Chapman and Hall, 1996).

Suwa, H.

H. Suwa and S. Todo, “Markov chain Monte Carlo method without detailed balance,” Phys. Rev. Lett. 105, 120603 (2010).
[Crossref]

Tarvainen, T.

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

Teller, A. H.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
[Crossref]

Teller, E.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
[Crossref]

Todo, S.

H. Suwa and S. Todo, “Markov chain Monte Carlo method without detailed balance,” Phys. Rev. Lett. 105, 120603 (2010).
[Crossref]

Vauhkonen, M.

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

Yamada, Y.

Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev. 21, 185–205 (2014).
[Crossref]

Yamashita, O.

Yodh, A. G.

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Zhao, H.

H. Zhao and R. J. Cooper, “Review of recent progress toward a fiberless, whole-scalp diffuse optical tomography system,” Neurophotonics 5, 011012 (2017).
[Crossref]

Zhu, Q.

Q. Zhu and S. Poplack, “A review of optical breast imaging: multi-modality systems for breast cancer diagnosis,” Eur. J. Radiol. 129, 109067 (2020).
[Crossref]

Zuev, K. M.

L. S. Katafygiotis and K. M. Zuev, “Geometric insight into the challenges of solving high-dimensional reliability problems,” Prob. Eng. Mech. 23, 208–218 (2008).
[Crossref]

Appl. Sci. (1)

Y. Jiang, Y. Hoshi, M. Machida, and G. Nakamura, “A hybrid inversion scheme combining Markov chain Monte Carlo and iterative methods for determining optical properties of random media,” Appl. Sci. 9, 3500 (2019).
[Crossref]

Eur. J. Radiol. (1)

Q. Zhu and S. Poplack, “A review of optical breast imaging: multi-modality systems for breast cancer diagnosis,” Eur. J. Radiol. 129, 109067 (2020).
[Crossref]

IEEE Trans. Geosci. Remote Sens. (1)

I. Langmore, A. B. Davis, and G. Bal, “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. 51, 2903–2919 (2013).
[Crossref]

Inverse Probl. (4)

S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inverse Probl. 22, 175–195 (2006).
[Crossref]

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).
[Crossref]

S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009).
[Crossref]

V. A. Markel and J. C. Schotland, “On the convergence of the Born series in optical tomography with diffuse light,” Inverse Probl. 23, 1445–1465 (2007).
[Crossref]

Inverse Probl. Imaging (1)

G. Bal, I. Langmore, and Y. Marzouk, “Bayesian inverse problems with Monte Carlo forward models,” Inverse Probl. Imaging 7, 81–105 (2013).
[Crossref]

J. Chem. Phys. (1)

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953).
[Crossref]

J. Comput. Appl. Math. (1)

T. Ooura and M. Mori, “The double exponential formula for oscillatory functions over the half infinite interval,” J. Comput. Appl. Math. 38, 353–360 (1991).
[Crossref]

J. Opt. Soc. Am. (1)

J. Opt. Soc. Am. A (1)

J. Phys. Soc. Jpn. (1)

K. Hukushima and K. Nemoto, “Exchange Monte Carlo method and application to spin glass simulations,” J. Phys. Soc. Jpn. 65, 1604–1608 (1996).
[Crossref]

Med. Phys. (1)

R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32, 1128–1139 (2005).
[Crossref]

Neurophotonics (1)

H. Zhao and R. J. Cooper, “Review of recent progress toward a fiberless, whole-scalp diffuse optical tomography system,” Neurophotonics 5, 011012 (2017).
[Crossref]

Opt. Express (1)

Opt. Rev. (1)

Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev. 21, 185–205 (2014).
[Crossref]

Phys. Med. Biol. (2)

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).
[Crossref]

M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999).
[Crossref]

Phys. Rev. E (2)

V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E 70, 056616 (2004).
[Crossref]

T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E 58, 5355–5363 (1998).
[Crossref]

Phys. Rev. Lett. (1)

H. Suwa and S. Todo, “Markov chain Monte Carlo method without detailed balance,” Phys. Rev. Lett. 105, 120603 (2010).
[Crossref]

Prob. Eng. Mech. (2)

S.-K. Au and J. L. Beck, “Estimation of small failure probabilities in high dimensions by subset simulation,” Prob. Eng. Mech. 16, 263–277 (2001).
[Crossref]

L. S. Katafygiotis and K. M. Zuev, “Geometric insight into the challenges of solving high-dimensional reliability problems,” Prob. Eng. Mech. 23, 208–218 (2008).
[Crossref]

Quart. Appl. Math. (1)

K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quart. Appl. Math. 2, 164–168 (1944).
[Crossref]

SIAM J. Appl. Math. (1)

D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963).
[Crossref]

Other (4)

W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials (Academic, 1979).

W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov Chain Monte Carlo in Practice (Chapman and Hall, 1996).

R. Fletcher, “A modified Marquardt subroutine for nonlinear least squares,” Report AERE-R 6799 (The Atomic Energy Research Establishment, Harwell, 1971).

D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University, 2000).

Data Availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)

Fig. 1.
Fig. 1. (top) The value $a$ is plotted against Monte Carlo steps for the simulated annealing in Section 7. (bottom) The reconstruction is done with the Levenberg–Marquardt algorithm.
Fig. 2.
Fig. 2. (left) One absorber at the depth 10 mm. (right) One absorber at a deeper position at $y = 15\;{\rm mm} $.
Fig. 3.
Fig. 3. Two absorbers with separation of 20 mm reconstructed by the proposed method.
Fig. 4.
Fig. 4. Two absorbers with separation of 20 mm. (left) The reconstruction by the truncated SVD with 52 singular values. (right) The reconstruction by the truncated SVD with 80 singular values.

Equations (58)

Equations on this page are rendered with MathJax. Learn more.

{(1ctD0Δ+μa(r))u(r,t)=f(r,t),rΩ,t>0,D0νu(r,t)+1ζu(r,t)=0,rΩ,t>0,u(r,0)=0,rΩ,
{D0Δu(r)+μa(r)u(r)=f(r),rΩ,D0νu(r)+1ζu(r)=0,rΩ.
f(r,t)=g0δ(rrs(p))δ(t)orf(r)=g0δ(rrs(p)),
μa(r)=μ¯a
μa(r)=μ¯a+δμa(r),rΩROI.
{D0Δ(u(r)u0(r))+μ¯a(u(r)u0(r))=δμa(r)u(r),rΩ,D0ν(u(r)u0(r))+1ζ(u(r)u0(r))=0,rΩ.
u(r)=u0(r)ΩG(r,r)δμa(r)u(r)dr.
vk+1(r)=ΩG(r,r)δμa(r)vk(r)dr,k=0,1,,
uR=u0ev1/u0,uR2=u0ev1/u0exp[v2u012(v1u0)2].
u(r)uR(r)=u0(r)(1ev1(r)/u0(r))ΩG(r,r)δμa(r)uR(r)drΩG(r,r)δμa(r)(u(r)uR(r))dr.
1ev1(r)/u0(r)v1(r)u0(r)=1u0(r)ΩG(r,r)δμa(r)u0(r)dr,
ΩG(r,r)δμa(r)uR(r)dr=ΩG(r,r)δμa(r)u0(r)ev1(r)/u0(r)drΩG(r,r)δμa(r)u0(r)dr.
uuRL(Ω)C1δμaL(Ω)u0L(Ω)+C1δμaL(Ω)uuRL(Ω),
uuRL(Ω)C1δμaL(Ω)u0L(Ω),
ϕ(p)=lnu0(rd(p))u(rd(p)).
ϕR(p)=v1(rd(p))u0(rd(p)),ϕR2(p)=v1(rd(p))u0(rd(p))+12(v1(rd(p))u0(rd(p)))2v2(rd(p))u0(rd(p)).
|eϕ(p)|=|uuRu0+uRu0|C1δμaL(Ω)u0L(Ω)u0+ev1/u0.
S(r)=0,±1,±2,,±M2,rΩROI.
δμa(r)=ηfa¯(x)δ(yy0),
fa¯(x)=[a¯3+3(1+tanhx210)a¯2](1tanhx2),
a=a(min)+a(max)a(min)M+1(S+M2+1)(a(min),a(max)].
δμa(ri)=δμa(max)(SiM+12),Si=0,±1,±2,,±M2
G(r,t;r,s)=eμ¯ac(ts)4πD0(ts)e(xx)24D0c(ts)[e(yy)24D0c(ts)+e(y+y)24D0c(ts)4πD0c(ts)e(y+y)24D0c(ts)e(y+y+2D0c(ts)/4D0c(ts))2×erfc(y+y+2D0c(ts)/4D0c(ts))]
u0(r,t)=g0eμ¯act2πD0te(xxsi)2+y24D0ct×[1πD0cte(y+2D0ct/4D0ct)2erfc(y+2D0ct/4D0ct)].
g(y,t;y,s)=14πD0(ts)e(y+y)24D0c(ts)[1+e(y+y)2(yy)24D0c(ts)4πD0c(ts)e(y+y2D0c(ts)+D0c(ts))2×erfc(y+y2D0c(ts)+D0c(ts))].
0tG(r,t;r,s)u0(r,s)ds=eμ¯act0te(xx)24D0c(ts)e(xxs(p))24D0csg(y,t;y,s)g(y,s;0,0)ds.
uR(r,t;rs(p))=u0(r,t;rs(p))exp[eμ¯actu0(r,t;rsi)0δμa(r)(0te(xx)24D0c(ts)e(xxs(p))24D0csg(y,t;y,s)g(y,s;0,0)ds)dxdy].
uR(rd(p),t;rs(p);a)=u(rd(p),t)=u0(rd(p),t;rs(p))exp[ηeμ¯actu0(rd(p),t;rs(p))0tg(0,t;y0,s)g(y0,s;0,0)(fa(x)e(xd(p)x)24D0c(ts)e(xxs(p))24D0csdx)ds],
u0=g0eμ¯act2πD0te(xd(p)xs(p))24D0ct[1πD0cte(D0ct)2erfc(D0ct)].
ϕR(p)=ηeμ¯actu0(rd(p),t;rs(p))0tg(0,t;y0,s)g(y0,s;0,0)(fa(x)e(xd(p)x)24D0c(ts)e(xxs(p))24D0csdx)ds.
H(S)=12p=1MSDj=1Mt|Φ(p)ϕ~R(p)|2+α|SS(0)|,
h(t)=1texp(y024Dct)[1πDcte(y02Dct+Dct)2erfc(y02Dct+Dct)].
u(rd(p),tj;rs(p);a¯)u(rd(p),tj;rs(p);a)=ηeμa0ctj(2πD)20tjh(tjs)h(s)[da¯(a,x)(1tanhx2)e(xd(p)x)24Dc(tjs)e(xxs(p))24Dcsdx]ds,
ξ(a,x)=a2[a+3(1+tanhx210)].
Kp,i=δμa(max)|ω|G(rd(p),ri)G(ri,rs(p))G(rd(p),rs(p)),
ϕR(p)ϕ~R(p),ϕR2(p)ϕ~R2(p),
ϕ~R(p)=i=1NKp,i(SiM+12),
ϕ~R2(p)=i=1NKp,i(SiM+12)+12i1=1Ni2=1NKp,i1Kp,i2(Si1M+12)(Si2M+12)|ω|2i1=1Ni2=1N[G(rd(p),rs(p))]1×G(rd(p),ri1)δμa(ri1)G(ri1,ri2)δμa(ri2)G(ri2,rs(p)).
Ψ(S)=12p=1MSD|Φ(p)ϕ~R2(p)|2+αi=1N|SiSi(0)|,
g1=Kp,i1Kp,i2,g2=g12(δμa(max))2|ω|2G(rd(p),ri1)G(ri1,ri2)G(ri2,rs(p))G(rd(p),rs(p)).
|ϕ(p)g2|=|g1||ϕ(p)g2g1|=|g1||ϕ(p)||12G(rd(p),rs(p))G(ri1,ri2)G(rd(p),ri2)G(ri1,rs(p))|.
|g1||ϕ(p)g2|
H(S)=i=1Nj=1NJijSiSji=1NhiSi,
Jij=12M2p=1MSDKp,iKp,j,hi=Mj=1NJij+1M(p=1MSDΦ(p)Kp,iα).
Z={Si}eβH(S),
π(S)=eβH(S)Z.
π(S)π(S)=eβ(H(S)H(S))=exp[β(Jii(Si2Si2)+(2j=1jiNJijSj+hi)(SiSi))].
T10int(log10T)2T.
a¯=1.5.
tj=jΔt(j=1,,Mt),Δt=5,Mt=500.
a0=0.01.
G(r,r)=12πD00cos(q(xx))λ×(eλ|yy|+λ1λ+1eλ(y+y))dq,
xs(p)=±2,±6,,±30mm,xd(p)=0,±4,±8,,±28mm.
μ¯a=0.02mm1,D0=0.33mm1.
δμa=0.2mm1.
ζ=21+rd1rd,
rd=1.4399n2+0.7099n1+0.6681+0.0636n.
α=0.01,Thigh=105,Tlow=1010,M=256.

Metrics