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

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$):

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

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

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]:

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

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

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

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

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

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

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

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

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

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

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

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

Let us define

Then, we have

Therefore, we obtain

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

We obtain

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

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

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

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

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})$:

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

We have

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$),

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

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

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

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:

In this paper, we decrease $T$ as

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

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

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

In Fig. 1, we compare the SA developed in this paper with the Levenberg–Marquardt algorithm [23–25], 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.

 

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 | PDF

 

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 | 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)$.

 

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

Download Full Size | PDF

 

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 | PDF

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

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

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

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]

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:

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.