Optical tomography with discretized path integral

Bingzhi Yuan; Toru Tamaki; Takahiro Kushida; Yasuhiro Mukaigawa; Hiroyuki Kubo; Bisser Raytchev; Kazufumi Kaneda

doi:10.1117/1.JMI.2.3.033501

13 August 2015 Optical tomography with discretized path integral

Bingzhi Yuan; Toru Tamaki; Takahiro Kushida; Yasuhiro Mukaigawa; Hiroyuki Kubo; Bisser Raytchev; Kazufumi Kaneda;

Author Affiliations +

J. of Medical Imaging, 2(3), 033501 (2015). doi:10.1117/1.JMI.2.3.033501

Abstract

We present a framework for optical tomography based on a path integral. Instead of directly solving the radiative transport equations, which have been widely used in optical tomography, we use a path integral that has been developed for rendering participating media based on the volume rendering equation in computer graphics. For a discretized two-dimensional layered grid, we develop an algorithm to estimate the extinction coefficients of each voxel with an interior point method. Numerical simulation results are shown to demonstrate that the proposed method works well.

Yuan, Tamaki, Kushida, Mukaigawa, Kubo, Raytchev, and Kaneda: Optical tomography with discretized path integral

1. Introduction

Optical tomography¹2.3.4.5.6.7.^–⁸ is known as a safer alternative to x-ray tomography. Usually, tomography consists of a light source generating penetrative light and a detector capturing the light, which allows to estimate the inside of the object through which the light is passing. The most important application is x-ray computed tomography (CT), where x rays are used due to their penetrative property. The balance between the radiation exposure of the human body and the quality of the obtained results has been debated since the early days when x-ray CT was invented. Therefore, there is an urgent demand for a safer medical tomography, such as optical tomography.

Modeling the behavior of light plays an important role in optical tomography, and in the mesoscale, in which the wavelength of light is close to the scale of tissue, the radiative transport equation (RTE) is used for describing the behavior of light scattering.⁵^,⁹ At the macroscale,⁶ the time-independent or dependent RTE is often approximated with a diffusion equation.

Similarly, the computer graphics community has the used time-independent RTE, and in contrast to the (surface) rendering equation,¹⁰^,¹¹ often calls it the volume rendering equation (VRE).¹⁰^,¹²

(1)

(ω \cdot \nabla) L (x, ω) = - σ_{t} (x) L (x, ω) + σ_{s} (x) \int_{S^{2}} f_{p} (x, ω, ω^{'}) L (x, ω^{'}) d ω^{'},

and the notations will be introduced in the following sections. The use of VRE enables us to render volumes of participating media, such as fog, cloud, and fire through which light is penetrating, and to obtain realistic volume rendering images of such scenes.¹³^,¹⁴ The path integral, which can be considered as a discrete version of the continuous Feynman path integral,¹⁵^,¹⁶ has been recently employed to solve the VRE in an efficient way with Monte Carlo integration, such as Metropolis light transport¹⁷^,¹⁸ or bidirectional path tracing.¹⁹

In this paper, we propose an optical tomography method using path integral as a forward model and solving a nonlinear inverse problem that minimizes the discrepancy between measurements and model predictions in a least-squares sense. To the best of our knowledge, the discretized path integral has not been used in optical tomography before. In our work, we simplify the path integral with some assumptions. The path integral, as the name suggests, gathers (or integrates) the contributions of all possible paths of light.¹⁷^,¹⁸^,²⁰21.22.^–²³ We approximate the integral of an infinite number of paths with the sum of a finite number of paths, discretize a continuous medium into voxels of a regular grid, and continuous light paths into discrete ones (i.e., polylines). We deal with anisotropic scattering having a peak in the forward direction, which is different from other discretization methods using discrete ordinate or spherical harmonics.¹³^,²⁴^,²⁵ In this work, we focus on estimating the spatially varying extinction coefficient $σ_{t} (x)$ at each discretized voxel location of the medium while fixing scattering properties (e.g., scattering coefficients $σ_{s}$ and phase functions $f_{p}$ ). By separating the scattering properties from our problem, we formulate optical tomography as an optimization problem with inequality constraints solved by an interior point method.

An interior point method²⁶ is an iterative method to solve an optimization problem with inequality constraints describing a feasible region in which the optimal solution must reside. To this end, a series of nonconstrained optimization problems are constructed by combining the constraints and the original objective function and are solved by an ordinal gradient-based (Quasi-Newton) method.

To summarize our contribution, we reformulate the problem of optical tomography by combining a path integral with several simplifying assumptions to model the light transport in the participating media. This paper is an extension of our previous conference version²⁷^,²⁸ with additional theoretical background and additional experiments and discussions, and is structured as follows. In Sec. 2, we briefly review previous work related to path integrals and optical tomography. In Sec. 3, we describe how to model the light transport in participating media and turn optical tomography into an optimization problem. In Sec. 4, we show how to solve the optimization problems. Section 5 reports some simulation results, and Sec. 6 concludes the paper.

2. Related Work

In this section, we briefly review related work on optical tomography and path integrals in computer graphics.

Optical tomography⁴^,⁵ (or inverse transport,⁶^,⁷ inverse scattering,²⁹ scattering tomography³⁰^,³¹) is a problem in medical imaging using light sources to reconstruct the optical properties of tissue from measurements (time-dependent or stationary, angular-dependent or independent) at the surface boundary. Analytically solving the RTE [Eq. (1)] with boundary conditions is difficult, however, and approximations, such as discrete ordinates and $N$ ’th-order spherical harmonics ( $P_{N}$ approximation), are often used and solved numerically by, for example, finite element methods (FEM) or finite difference methods. The famous diffuse approximation⁵^,⁶ (DA) is a $P_{1}$ (thus first-order) approximation with the assumption on a phase function being isotropic. The DA is an approximation to RTE at a macroscopic scale when scattering is large while absorption is low and scattering is not highly peaked. Diffuse optical tomography (DOT) is based on DA and today represents the frontier of optical tomography³²^,³³ with many clinical applications.³⁴ DA, however, does not often hold in realistic participating (scattering) media; absorption may not be small compared to scattering, and the shapes of the phase functions can be highly peaked in the forward direction, which is often modeled by Henyey-Greenstein,³⁵ Schlick,³⁶ or Mei and Rayleigh phase functions.¹⁰^,¹²^,³⁷^,³⁸ Experimental evidence³⁹ also suggests a highly peaked shape of the phase functions in biological media. DOT works, but is still limited; therefore, other methods have also been studied for cases when DA does not hold.

Statistical Monte Carlo methods are used for media in which the assumptions do not hold; however, they are computationally intensive and inefficient for solving the forward problem,⁴5.6.^–⁷^,³⁴ i.e., solving the RTE with given parameters. Therefore, Monte Carlo based approaches have been used for estimating the spatially constant (not varying) parameters in homogeneous media, such as paper,⁴⁰^,⁴¹ clouds,⁴² liquids,⁴³ plastics,⁴⁴ or uniform material samples.⁴⁵ Another difficulty of Monte Carlo based inverse methods is that an analytical forward model prediction is hard to obtain when we want to minimize the difference between the prediction and measurements except for very special structures.⁴⁶^,⁴⁷ A gradient based least-square approach has been proposed but only for spatially constant parameter estimation,⁴⁰^,⁴¹^,⁴⁸ while model-free approaches have relied on genetic algorithms,⁴²^,⁴⁴ numerical perturbation,⁴⁹^,⁵⁰ voting,⁵¹ or even simple backprojection.⁵² One of the contributions of the current paper is to enable us to use a gradient based optimization approach for estimating spatially varying parameters, which is extensible by using many optimization methods.

Similar to optical tomography, modeling light transport plays a very important role in computer graphics. Our own work on optical tomography is inspired by Monte Carlo based statistical methods. In the last two decades, methods based on path integrals¹⁷18.^–¹⁹^,⁵³54.^–⁵⁵ have provided models of light transport for efficient volume rendering. For solving RTE, a path integral has been used for a forward problem solver,¹⁶^,⁵⁶^,⁵⁷ and has also been applied to optical tomography, but under the diffusion assumption.⁵⁸^,⁵⁹ Our proposed method is based on a path integral to explicitly express the forward model prediction, which is very suitable for solving the inverse problem with gradient based methods. This is an advantage of our method over existing methods because the paths used in the forward model can be generated by either a deterministic or statistical (Monte Carlo) method. To achieve an efficient forward model, we introduce a simplified layered scattering model that uses a limited number of deterministic paths instead of Monte Carlo simulated ones.

3. Method: Forward Problem

We deal with the following optical tomography problem [this is a conceptual formulation and the actual problem is shown in Eq. (29)].

(2)

\min_{σ_{t}} \sum_{i, j} {| I_{i j} - P_{i j} (σ_{t}) |}^{2},

where

σ_{t}

is a vector representing the spatial distribution of the extinction coefficients to be estimated. We divide our discussion into two parts: forward and inverse problems. The forward problem focuses on building a mathematical model

P_{i j} (σ_{t})

of the light transport between a light source

i

and a detector

j

. We will make some assumptions on the light transport and the medium to simplify the forward model. An inverse problem minimizes the difference between the observations

I_{i j}

of the detector and the forward model to estimate the spatial distribution of the extinction coefficients

σ_{t}

.

3.1.

Forward Model

In the forward problem, as we mentioned before, we use a path integral to build a mathematical model for the light transport. Here, we follow the notation developed in the computer graphics literature¹⁷^,²³^,⁵³^,⁶⁰ to introduce the path integral. Sections 3.2 to 3.6 will show the simplified model we propose.

Given a space $R^{3}$ , a light source is located at $x_{0} \in R^{3}$ and a detector at $x_{M + 1} \in R^{3}$ , and in between them is the participating media $ν \subset R^{3}$ with boundary $\partial ν$ and interior volume $ν_{0} ≔ ν ∖ \partial ν$ . A light path $\tilde{x}$ connecting $x_{0}$ and $x_{M + 1}$ of length $M + 2$ consists of $M + 2$ vertices $x_{m} \in R^{3}$ for $m = 0, 1, \dots, M + 1$ , denoted by $\tilde{x} = x_{0}, x_{1}, \dots, x_{M}, x_{M + 1}$ . Thus, absorption, scattering, or reflection events happen at $x_{1}, \dots, x_{M}$ . The set of all paths of length $M$ is denoted by $Ω_{M}$ . The path space $Ω$ is the countable set of all paths $Ω_{M}$ of finite length.

(3)

Ω = ⋃_{M = 2}^{\infty} Ω_{M} .

A direction is denoted by $ω \in S^{2}$ , where $S^{2}$ is a unit sphere in $R^{3}$ . A unit vector $ω_{x_{m}, x_{m + 1}}$ is the direction from vertex $x_{m}$ to vertex $x_{m + 1}$ in a path $\tilde{x}$ .

Veach²⁰ introduced a framework representing the rendering equation in the form of a path integral for scenes without participating media (i.e., no scattering), and later, Pauly et al.¹⁷ extended it to the volume rendering equation with scattering. The amount of light $I$ observed by the detector is given by the path integral

(4)

I = \int_{Ω} f (\tilde{x}) d μ (\tilde{x}),

which is an integral over the path space. Here,

μ (\tilde{x})

is a measure of path

\tilde{x}

.

(5)

d μ (\tilde{x}) = \prod_{m = 0}^{M + 1} d μ (x_{m}), d μ (x_{m}) = {\begin{cases} d A (x_{m}), & x_{m} \in \partial ν \\ d V (x_{m}), & x_{m} \in ν_{0} \end{cases},

where

d μ (x_{m})

denotes the differential measure at vertex

x_{m}

.

f (\tilde{x})

is a measurement contribution function defined as follows:

(6)

f (\tilde{x}) = L_{e} (x_{0}, x_{1}) G (x_{0}, x_{1}) [\prod_{m = 1}^{M} f_{f} (x_{m - 1}, x_{m}, x_{m + 1}) G (x_{m}, x_{m + 1})] W_{e} (x_{M}, x_{M + 1}),

where

W_{e} (x_{M}, x_{M + 1})

is the camera response function, and

L_{e} (x_{0}, x_{1})

is the intensity of the light emitted from the light source

x_{0}

to vertex

x_{1}

.

f_{f} (x_{m - 1}, x_{m}, x_{m + 1})

is a scattering kernel at

x_{m}

with respect to the locations of vertices

x_{m - 1}

and

x_{m + 1}

.

(7)

f_{f} (x_{m - 1}, x_{m}, x_{m + 1}) = {\begin{cases} f_{s} (x_{m - 1}, x_{m}, x_{m + 1}), & x_{m} \in \partial ν \\ σ_{s} (x_{m}) f_{p} (x_{m - 1}, x_{m}, x_{m + 1}), & x_{m} \in ν_{0} \end{cases} .

Here, the bidirectional scattering distribution function $f_{s} (x_{m - 1}, x_{m}, x_{m + 1})$ is used for locations on the surface of objects, and the scattering coefficient $σ_{s} (x_{m})$ at $x_{m}$ and the phase function $f_{p} (x_{m - 1}, x_{m}, x_{m + 1})$ are used for those inside the medium. $G (x_{m}, x_{m + 1})$ is a generalized geometric term.

(8)

G (x_{m}, x_{m + 1}) = T (x_{m}, x_{m + 1}) g (x_{m}, x_{m + 1}),

where

g (x_{m}, x_{m + 1})

is a geometric term.

(9)

g (x_{m}, x_{m + 1}) = {\begin{cases} \frac{| n_{g} (x_{m}) \cdot ω_{x_{m}, x_{m + 1}} |}{{‖ x_{m} - x_{m + 1} ‖}^{2}}, & x_{m} \in \partial ν \\ \frac{1}{{‖ x_{m} - x_{m + 1} ‖}^{2}}, & x_{m} \in ν_{0} \end{cases},

with unit normal

n_{g} (x_{m})

of the surface at

x_{m} \in \partial ν

.

T (x_{m}, x_{m + 1})

is a transmittance that describes the attenuation when light passes through the medium.

(10)

T (x_{m}, x_{m + 1}) = {\begin{cases} e^{- τ (x_{m}, x_{m + 1})}, & {x_{m}, x_{m + 1}} \subset ν_{0} \cup \partial ν \\ 0, & otherwise \end{cases},

(11)

τ (x_{m}, x_{m + 1}) = \int_{0}^{1} σ_{t} [(1 - s) x_{m} + s x_{m + 1}] d s,

where

σ_{t} (x_{m})

is the extinction coefficient at vertex

x_{m}

.

Putting all together, we have a path integral of the following infinite sum of all possible path contributions.

(12)

I = \sum_{M = 2}^{\infty} \sum_{k \in Ω_{M}} L_{e} (x_{0}, x_{1}) G (x_{0}, x_{1}) [\prod_{m = 1}^{M} f_{f} (x_{m - 1}, x_{m}, x_{m + 1}) G (x_{m}, x_{m + 1})] W_{e} (x_{M}, x_{M + 1}) \prod_{m = 0}^{M + 1} d μ (x_{m}) .

Note that all vertices ${x_{m}}$ depend on a path $k$ ; different paths have different sets of vertices. In the equation above, however, we omit the path index $k$ for simplicity. Later, we will again use $k$ as the path index.

3.2.

Assumptions on the Path Integral Formulation

As our target is optical tomography, we restrict the model to deal with inside the participating media. To do so, we assume that the light source $x_{0}$ and detector $x_{M + 1}$ are located on the surface, and the other vertices $x_{1}, x_{2}, \dots, x_{M}, x_{M + 1}$ are inside the medium, that is, $x_{0}, x_{M + 1} \in \partial ν$ and $x_{1}, \dots, x_{M} \in ν_{0}$ . Then the transmittance is simplified as

(13)

T (x_{m}, x_{m + 1}) = e^{- τ (x_{m}, x_{m + 1})} .

Furthermore, we assume that the observations are ideal and the camera response function is the identity, $W_{e} (x_{M}, x_{M + 1}) = 1$ .

Apart from the assumptions above, we rewrite the geometric term and the differential measure. The definitions above use area measures $d A (x_{m})$ and volume measures $d V (x_{m})$ along with the squared distance geometric term;¹⁷^,²³^,⁵³ however, steradian measures $d ω (x_{m})$ and the identity geometric term is equivalent and also widely used.¹⁰^,¹²^,⁶⁰

(14)

g (x_{m}, x_{m + 1}) d μ (x_{m}) = d ω (x_{m}) .

Therefore, we employ the steradian measures and rewrite it as follows:

(15)

g (x_{m}, x_{m + 1}) = 1,

(16)

d μ_{k} (x_{m}) = {\begin{cases} d A (x_{0}), & m = 0 \\ d ω (x_{m}), & m = 1, \dots, M + 1 \end{cases} .

Now, Eq. (12) is written as

(17)

I = \sum_{M = 2}^{\infty} \sum_{k \in Ω_{M}} L_{e} (x_{0}, x_{1}) T (x_{0}, x_{1}) d A (x_{0}) [\prod_{m = 1}^{M} f_{f} (x_{m - 1}, x_{m}, x_{m + 1}) T (x_{m}, x_{m + 1}) d ω (x_{m})] d ω (x_{M + 1}) .

3.3.

Discretization of the Forward Model

For numerical computation, we first discretize the medium into voxels of a regular grid, where each voxel has its own extinction coefficient $σ_{t} [b]$ ( $b$ is the index of the voxel) as shown in Fig. 1.

Fig. 1

Illustration of a discretization example. (a) Voxelization of the medium into a regular grid of size $5 \times 5$ . Voxels are indexed in raster scan order in this example, from left to right, and top to bottom. Each voxel $b$ has extinction coefficient $σ_{t} [b]$ . (b) A path segment between vertices $x_{1}$ and $x_{2}$ . Voxels involved in the segment are shaded. (c) Lengths $d_{12} [b]$ of the involved voxels $b = 2, 3, 8, 9$ . Here we denote $d_{12} [b]$ instead of $d_{x_{1}, x_{2}} [b]$ for simplicity.

With this voxelization, the paths of light are also divided into segments, as explained below. First, we explain the integral [Eq. (11)] along a single segment $x_{m}, x_{m + 1}$ of a path $\tilde{x}$ . It describes the attenuation of light along the segment due to the extinction coefficients of the voxels involved. Because of the discretization of the medium, Eq. (11) can be written as a sum of voxel-wise multiplications.

(18)

τ (x_{m}, x_{m + 1}) = \int_{0}^{1} σ_{t} [(1 - s) x_{m} + s x_{m + 1}] d s = \sum_{b \in B_{x_{m}, x_{m + 1}}} σ_{t} [b] d_{x_{m}, x_{m + 1}} [b] = σ_{t}^{T} d_{x_{m}, x_{m + 1}} .

For the second equality, $b$ is the index of a set $B_{x_{m}, x_{m + 1}}$ of all voxels involved by segment $x_{m}, x_{m + 1}$ , and $d_{x_{m}, x_{m + 1}} [b]$ is the length of the part of the segment $x_{m} x_{m + 1}$ passing through voxel $b$ . This is illustrated in Fig. 1(c). The extinction coefficient $σ_{t}$ is now a piece-wise constant function because of the voxelization; then the integral turns into a sum (the idea that this integral can be turned into a sum has been discussed before,⁶¹ however, not in the context of tomography).

This simplifies the computation; however, the sum over a set $B_{x_{m}, x_{m + 1}}$ is not preferable in terms of implementation and optimization. We propose here to use a vector representation of both extinction coefficients and segment lengths, which is the third equality of the above equation. The first vector $σ_{t}$ stores the values of the extinction coefficients $σ_{t} [b]$ of all voxels. This vector can be generated by serializing the voxels on the grid in a certain order. The second vector $d_{x_{m}, x_{m + 1}}$ contains the values of the lengths $d_{x_{m}, x_{m + 1}} [b]$ for all voxels. We should note that this vector is very sparse; most of the voxels have no intersection with the segment $x_{m}, x_{m + 1}$ . Hence, only a few elements in $d_{x_{m}, x_{m + 1}}$ have nonzero values, and the other elements are zero because those voxels $b$ have no intersection and $d_{x_{m}, x_{m + 1}} [b] = 0$ .

This sparsity of the vector facilitates the construction of a whole path $\tilde{x}$ because path segments can be added as follows:

(19)

D_{k} = \sum_{m = 0}^{M} d_{x_{m}, x_{m + 1}},

where

D_{k}

is the vector of a complete path

k

of length

M + 2

; the

b

’th element can be interpreted as the length of the segment when the path passes through voxel

b

. This notation simplifies a part of Eq. (17) as follows:

(20)

\prod_{m = 0}^{M} T (x_{m}, x_{m + 1}) = \prod_{m = 0}^{M} e^{- τ (x_{m}, x_{m + 1})} = e^{- \sum_{m = 0}^{M} τ (x_{m}, x_{m + 1})} = e^{- \sum_{m = 0}^{M} σ_{t}^{T} d_{x_{m}, x_{m + 1}}} = e^{- σ_{t}^{T} D_{k}} .

Using this notation to rewrite Eq. (17), we have

(21)

I = \sum_{M = 2}^{\infty} L_{e} (x_{0}, x_{1}) \sum_{k \in Ω_{M}} H_{k} e^{- σ_{t}^{T} D_{k}} = L_{e} (x_{0}, x_{1}) \sum_{k \in Ω} H_{k} e^{- σ_{t}^{T} D_{k}},

where the factor

H_{k}

, defined as

(22)

H_{k} = d A (x_{0}) d ω (x_{M + 1}) \prod_{m = 1}^{M} f_{f} (x_{m - 1}, x_{m}, x_{m + 1}) d ω (x_{m}),

describes the contributions of the scattering coefficients and phase functions, and the exponential factor represents attenuation due to absorption (and outscattering) over the path.

3.4.

Two-Dimensional Layered Model of Forward Scattering

As a first attempt, we design a two-dimensional (2-D) layered grid, instead of the three-dimensional (3-D) one. Since we voxelize the medium into a regular grid, the 2-D medium consists of parallel layers. Hereafter, a 3-D direction $ω$ between vertices is written as a 2-D direction $θ$ and a steradian measure $d ω$ as an angular measure $d θ$ .

As shown in Fig. 2, we assume a particular layer scattering having the following properties. First, vertices $x_{1}, \dots, x_{M}$ of path $\tilde{x}$ are located at the centers of each voxel. The light source $x_{0}$ is located on the boundary of the top surface of the voxels in the top layer. Similarly, the detector $x_{M + 1}$ is located on the boundary of the bottom surface of the voxels in the bottom layer. Second, directions $θ_{x_{0}, x_{1}}$ and $θ_{x_{M}, x_{M + 1}}$ at the beginning and end of a path are perpendicular to the boundary. This means that scattering begins at $x_{1}$ and ends at $x_{M}$ . Third, forward scattering happens layer by layer. More specifically, light is scattered at the center of a voxel in a layer and then goes to the center of a voxel in the next (below) layer. Scattering is assumed to happen every time the light traverses voxel centers. Even if the next voxel is just below the current voxel and the path segment is straight, it is regarded as scattering. Fourth, the scattering coefficient is uniform, $σ_{s} (x) = σ_{s}$ .

Fig. 2

Proposed two-dimensional layered model of scattering. This example shows path $\tilde{x}$ consisting of vertices $x_{1}, \dots, x_{M}$ located at the centers of voxels in a grid with $M$ parallel layers. $x_{0}$ is a light source located on the top surface, and $x_{M + 1}$ is a detector at the bottom. At each vertex, the light scatters to voxels in the next layer, and possible scattering directions are indicated by arrows.

By ignoring paths exiting from the sides of the grid, the number of all possible paths is $N^{M}$ , where $M$ is the number of layers and $N$ is the number of voxels in one layer.

3.5.

Approximating the Phase Function with a Gaussian

We use a Gaussian model $f_{p} (θ, σ^{2})$ as an approximation of the phase function

(23)

f_{p} (x_{m - 1}, x_{m}, x_{m + 1}) \equiv f_{p} (θ_{m}, σ^{2}) = \frac{1}{\sqrt{2 π σ^{2}}} \exp (\frac{- θ_{m}^{2}}{σ^{2}}), - \frac{π}{2} < θ_{m} < \frac{π}{2},

where the variance

σ^{2}

controls the scattering property; larger values of

σ^{2}

mean strong forward scattering. This Gaussian approximation is convenient in our case because of the following two reasons.

First, existing phase function models¹⁰^,¹²^,³⁵36.37.^–³⁸ are those for 3-D scattering, not for 2-D. This means that those functions are normalized for integrals over the unit sphere $S^{2}$ : $\int_{S^{2}} f_{p} (ω) d ω = 1$ . Most of the phase functions assume isotropy (rotational symmetry), and hence, the function has a form taking angle $θ$ as an argument; however, $\int_{- π}^{π} f_{p} (θ) d θ \neq 1$ . These functions, therefore, are not adequate for our case.

Second, our assumption of layer-wise forward scattering does not allow scattering to happen backwards or sideways, and the Gaussian model is suitable for it. As shown in Fig. 3, the Gaussian model has the form of forward-only scattering (no backwards or sideways) in a reasonable range of $σ^{2}$ , and it is almost normalized; $\int_{- π / 2}^{π / 2} f_{p} (θ, σ^{2}) d θ \approx 1$ . Other 2-D phase functions exist which are not forward-only. For example, Heino et al.⁶² introduced a 2-D analog of Henyey-Greenstein’s phase function,³⁵ shown in Fig. 3. Although the parameters are different, the two functions in Fig. 3 have similar shapes. The most important difference is that Heino’s function has backward scattering, but our Gaussian model does not. More realistic scattering rather than the layer-wise forward scattering introduced here needs Heino’s or Henyey-Greenstein’s phase function.

Fig. 3

Comparison of two-dimensional phase functions. The upward vertical direction is $θ = 0$ , and horizontal directions are $θ = \pm (π / 2)$ . (a) Gaussian approximated phase functions with $σ^{2} = 0.1, 0.2, \dots, 1.0$ . The tallest and narrowest shapes correspond to $σ^{2} = 0.1$ , and the shape becomes shorter and rounder for larger values of $σ^{2}$ . (b) Heino’s two-dimensional analogs⁶² of Henyey-Greenstein’s phase function with parameter $g = 0.1, 0.2, \dots, 1.0$ . The tallest and narrowest shapes correspond to $g = 1.0$ , and the shape becomes shorter and close to a hemisphere for smaller values of $g$ .

We should note one further simplification in our layer-wise forward scattering model. The angle $θ_{m}$ in the phase function is usually defined between $θ_{x_{m - 1}, x_{m}}$ and $θ_{x_{m}, x_{m + 1}}$ , that is, the difference of directions changed by the scattering event. Instead of dealing with such an exact difference of directions, we use the angle between $θ_{x_{m}, x_{m + 1}}$ and the vertical (downward) direction for efficiency of computation. This assumption enables us to discretize the Gaussian phase function much more easily. Since $f_{p} (θ)$ integrates to (approximately) one, such a normalization can be discretized with a sum as follows:

(24)

\int_{- π / 2}^{π / 2} f_{p} (θ, σ^{2}) d θ \approx \sum_{b \in B_{n}} f_{p} (θ_{b}, σ^{2}) Δ θ_{b} \approx 1,

where

B

is a set of voxel indices in the next layer

n

,

θ_{b}

is an alternative form of the corresponding

θ_{x_{m}, x_{m + 1}}

, and

Δ θ_{b}

is the angle measure as shown in Fig. 4.

Fig. 4

An illustration of angle measure $Δ θ_{b}$ for voxel $b$ in the next layer. For the center voxel of the upper layer, voxel $b$ (shaded) in the next layer subtends an angle of $Δ θ_{b}$ , which is used for the angle measure in Eq. (24).

The above equation can be considered as the energy distribution from a voxel in one layer to the voxels in the next layer. For a voxel $b$ at direction $θ_{b}$ , the value of $f_{p} (θ_{b}, σ^{2}) Δ θ_{b}$ describes what percentage of the energy will be scattered to this voxel. Figure 5 shows plots of the values corresponding to two phase functions with different parameters. We can see that, due to forward scattering, most of the energy is concentrated in the voxel just below, and a small part goes to the adjacent voxels.

Fig. 5

(a) The phase functions with parameter $σ^{2} = 0.2$ (dashed line) and $σ^{2} = 0.4$ (solid line). Plot of the value $f_{p} (θ_{b}, σ^{2}) Δ θ_{b}$ for each voxel $b$ for (b) $σ^{2} = 0.2$ and (c) $σ^{2} = 0.4$ . Note that index $b$ is relative to the voxel in the next layer just below the voxel in consideration. The voxel just below is $b = 0$ , the voxel on its right side is $b = 1$ , and that on the left side is $b = - 1$ .

The contribution $H_{k}$ in Eq. (22) now needs to be rewritten so that it deals with the Gaussian phase function and the discretized energy distribution discussed above. First, we reorder the measure

(25)

H_{k} = d A (x_{0}) d θ (x_{M + 1}) \prod_{m = 1}^{M} f_{f} (x_{m - 1}, x_{m}, x_{m + 1}) d θ (x_{m})

(26)

= d A (x_{0}) d θ (x_{1}) \prod_{m = 1}^{M} f_{f} (x_{m - 1}, x_{m}, x_{m + 1}) d θ (x_{m + 1}),

and then replace the factors with the Gaussian phase function.

(27)

H_{k} = d A (x_{0}) Δ θ_{x_{0}, x_{1}} σ_{s}^{M} \prod_{m = 1}^{M} f_{p} (θ_{x_{m}, x_{m + 1}}, σ^{2}) Δ θ_{x_{m}, x_{m + 1}} .

Note that the factor $d A (x_{0}) Δ θ_{x_{0}, x_{1}} σ_{s}^{M}$ is common for all paths because we assumed that the grid is uniform so that $d A (x_{0})$ is constant, and the direction $θ_{x_{0}, x_{1}}$ (or $ω_{x_{0}, x_{1}}$ ) is perpendicular to the top surface, and $σ_{s}$ is constant.

3.6.

Observation Model

Suppose the 2-D layered medium is an $M \times N$ grid; it has $M$ layers, each of which is made of $N$ voxels. We now construct an observation model of the light transport between a light source and a detector: emitting light to each of the voxels at the top layer, and capturing light from each voxel from the bottom layer. More specifically, let $i \in B_{1}$ and $j \in B_{M}$ be voxel indices of the light source and detector locations, respectively. By restricting the light paths to only those connecting $i$ and $j$ , the observed light $I_{i j}$ is written as follows:

(28)

I_{i j} = I_{0} \sum_{k = 1}^{N_{i j}} H_{i j k} e^{- σ_{t}^{T} D_{i j k}},

where

H_{i j k}

and

D_{i j k}

are the same as in Eqs. (27) and (21), respectively, but are restricted to paths connecting

i

and

j

, and

I_{0} = L_{e} (x_{0}, x_{1})

assuming the light source to be constant.

In the above equation, $k$ indexes the light paths, which share the same $i$ and $j$ . Due to the layered scattering model in the $N \times M$ grid, the number of different paths between $i$ and $j$ is $N_{i j} = N^{M - 2}$ . This is, however, too large even for small $N$ and $M$ , e.g., $N = M = 10$ . Therefore, we exclude paths having small contributions from the computation. This is done by a simple thresholding while computing $H_{i j k}$ as shown in Algorithm 1. This results in generating fewer paths; $N_{i j} \leq N^{M - 2}$ . For example, there are $N_{i j} = 742$ paths for $N = M = 20$ with $σ^{2} = 0.4$ when $th = 0.001$ , which enables us to reduce the computation cost.

Algorithm 1

Computing contribution Hijk and omitting low contribution path by thresholding.

Input: Threshold th, path x˜=x0,⋯,xM+1.

Output: Contribution Hijk.

1 Hijk=1;

2 form=1toMdo

3 Hijk=Hijkfp(θxm,xm+1,σ2)Δθxm,xm+1

4 ifHijk≤ththen

5 stop;

6 omit this path;

7 accept this path;

8 return Hijk;

4. Method: Inverse Problem

Next, we propose a method for the inverse problem of the forward model [Eq. (28)] to estimate the extinction coefficients of the 2-D layered model. As we mentioned before, we fix the light paths and assume that the scattering coefficients and parameters of the Gaussian phase function are uniform and known in advance.

4.1.

Cost Function

In the $M \times N$ 2-D layered medium described in Sec. 3.6, we had assumed a configuration of a light source and detector similar to the left-most one shown in Fig. 6; the light source is located above the medium and the detector is below, and the observed light is $I_{i j}$ , where $i, j$ are the voxel indices of the light source and detector locations. By sliding the light source and the detector, we can obtain $N^{2}$ observations, resulting in the following least-squares equation:

(29)

\min_{σ_{t}} f_{0}, f_{0} = \sum_{i = 1}^{N} \sum_{j = 1}^{N} {| I_{i j} - I_{0} \sum_{k = 1}^{N_{i j}} H_{i j k} e^{- σ_{t}^{T} D_{i j k}} |}^{2},

under

2 M N

constraints

(30)

0 ⪯ σ_{t} ⪯ u,

where

⪯

denotes the generalized inequality, i.e., all elements in the vector must satisfy the inequality. The lower bound 0 comes from the fact that any media must have positive extinction coefficients, while the upper bound

u

is used for numerical stability to exclude unrealistic values to be estimated.

Fig. 6

Four configurations of light sources and detectors. From left to right, we call configurations T2B (top-to-bottom), L2R (left-to-right), B2T (bottom-to-top), and R2L (right-to-left), which represent locations of light sources and detectors.

Furthermore, as shown in Fig. 6, we have four configurations of light sources and detectors by changing their positions. This gives us four different sets of observations $I_{i j}$ and paths $i j k$ . These four different sets lead to four objective functions ( $f_{T 2 B}$ , $f_{L 2 R}$ , $f_{B 2 T}$ , $f_{R 2 L}$ ) as shown in Fig. 6. Since the four objective functions share the same variables $σ_{t}$ , we can use all of them at the same time by adding them to form a new single function $f_{0}$ at the expense of additional (factor of four) computation cost.

(31)

\min_{σ_{t}} f_{0}, f_{0} = f_{0}^{T 2 B} + f_{0}^{L 2 R} + f_{0}^{B 2 T} + f_{0}^{R 2 L} subject to 0 ⪯ σ_{t} ⪯ u .

4.2.

Optimization Problem with Inequality Constraints

Since the inverse problem [Eq. (31)] is nonlinear, we employ an interior point method,²⁶ an iterative optimization algorithm for problems with constraints. Here, we first review several key points in optimization; then we will develop an algorithm to solve Eq. (31) along with the required first- and second-order derivatives of the cost function.

4.2.1.

Unconstrained problem: Quasi-Newton

First, we review optimization without constraints, which is used inside the interior point method. The general form of unconstrained optimization is

(32)

\min_{σ_{t}} f (σ_{t}),

where

σ_{t} \in R^{N \times M}

is a real vector and

f : R^{N \times M} \to R

is an objective function that is twice continuously differentiable.

To solve it, an iterative procedure begins with an initial guess $σ_{t}^{0}$ and generates a sequence ${σ_{t}^{k}}_{k = 0}^{\infty}$ . It stops when the change of solutions is small enough. The information about function $f$ at $σ_{t}^{k}$ or even previous estimates $σ_{t}^{0}, σ_{t}^{1}, \dots, σ_{t}^{k - 1}$ is used to calculate a direction $p_{k}$ to move with a step size $α_{k}$ . A line search is often used to determine the step size by searching along the direction starting from $σ_{t}^{k}$ for finding $σ_{t}^{k + 1}$ with the least value of the objective function

(33)

\min_{α_{k} > 0} f (σ_{t}^{k} + α_{k} p_{k}) .

Once we find the step size, the estimate $σ_{t}^{k + 1}$ is updated as $σ_{t}^{k + 1} \leftarrow σ_{t}^{k} + α_{k} p_{k}$ . The direction is $p_{k} = - B_{k} \nabla f (σ_{t}^{k})$ for the Newton’s method, where $B_{k} = \nabla^{2} f {(σ_{t}^{k})}^{- 1}$ is the inverse of the Hessian.

The Newton’s method is well known for its second-order convergence and accuracy. However, when the dimension of the problem is large, calculating the Hessian and its inverse is computationally expensive. Therefore, Quasi-Newton methods are often used, where the inverse Hessian is updated by incremental approximations in order to reduce the computation cost. The Broyden-Fletcher-Goldfarb-Shanno (BFGS) update rules are well known.⁶³

(34)

s = σ_{t}^{k} - σ_{t}^{k - 1},

(35)

y = \nabla f (σ_{t}^{k}) - \nabla f (σ_{t}^{k - 1}),

(36)

B_{k} = (I - \frac{s y^{T}}{y^{T} s}) B_{k - 1} (I - \frac{y s^{T}}{y^{T} s}) + \frac{s s^{T}}{y^{T} s} .

When the conditions $y^{T} s > 0$ and $B_{0} ≻ 0$ (where $≻ 0$ means positive definite) are satisfied, the BFGS update guarantees the positive definiteness of $B_{k}$ . Algorithm 2 shows the Quasi-Newton method.

Algorithm 2

The Quasi-Newton method with BFGS update rule.

Input: A feasible initial solution σt0, and B0≻0.

Result: An estimate σt⋆.

1 repeat

2 Compute the Quasi-Newton direction: pk=−Bk∇f(σtk).

3 Find step length αk with line search.

4 Update estimate σtk+1←σtk+αkpk.

5 Update Bk with BFGS.

6 until convergence;

4.2.2.

Constrained problem: interior point

Here we introduce a constrained optimization with inequality constraints of the form

(37)

\min_{σ_{t}} f_{0} (σ_{t}) subject to f_{i} (x) \leq 0, i = 1, \dots, m,

where

σ_{t} \in R^{N \times M}

is a real vector and

f_{0}, \dots, f_{m} : R^{N \times M} \to R

are twice continuously differentiable.

The idea is to approximate it as an unconstrained problem. Using Lagrange multipliers, we can first rewrite Eq. (37) as

(38)

\min_{σ_{t}} f_{0} (σ_{t}) + \sum_{i = 1}^{m} I [f_{i} (σ_{t})],

where

I : R \to R

is an indicator function, which keeps the solution inside the feasible region.

(39)

I (f) = {\begin{cases} 0, & f \leq 0 \\ \infty, & f > 0 \end{cases} .

Equation (38) now has no inequality constraints, while it is not differentiable due to $I$ .

The barrier method²⁶ is an interior point method that introduces a logarithmic barrier function to approximate the indicator function $I$ as follows:

(40)

\hat{I} (f) = - (1 / t) \log (- f),

where

t > 0

is a parameter to adjust the accuracy of approximation. The log barrier function goes to infinity rapidly as

f

goes close to 0, while it is close to 0 when

f

is far away from 0. Since

\hat{I} (f)

is differentiable, we have

(41)

\min_{σ_{t}} f_{0} (σ_{t}) + \sum_{i = 1}^{m} - (1 / t) \log [- f_{i} (σ_{t})],

or equivalently,

(42)

\min_{σ_{t}} t f_{0} (σ_{t}) - \sum_{i = 1}^{m} \log [- f_{i} (σ_{t})] .

The barrier method solves Eq. (42) iteratively by increasing the parameter $t$ . At the limit of $t \to \infty$ , the above problem coincides with the original problem [Eq. (38)].

4.3.

Algorithm for Solving the Inverse Problem

Algorithm 3 shows our algorithm, which uses a barrier method with Quasi-Newton for solving the inverse problem. We should mention the following parts where we have modified the original algorithm.²⁶

Algorithm 3

Barrier method of interior point with Quasi-Newton solver.

Data: Parameters μ>1, ϵ>0, and t=tinit>0.

Input: A feasible initial estimate σt0, and B≻0.

Result: An estimate σt⋆.

1 while(2MN/t)≥ϵdo // outer loop: barrier method

2 t←μt.

3 Set a log-barriered cost function

f(t)=tf0−∑b{log(σt[b])+log(u−σt[b])}

4 k←0, Bk←B, σtk←σt.

5 repeat // inner loop: Quasi-Newton

6 Compute the Quasi-Newton direction: pk=−Bk∇f(σtk).

7 Find step length αk with line search.

8 whileσtk+αkpk is not feasible do

9 Halve the step size: αk←αk/2.

10 Update estimate σtk+1←σtk+αkpk.

11 s=σtk+1−σtk.

12 y=∇f(σtk+1)−∇f(σtk).

13 ifyTs>0then

14 Update Bk+1 with BFGS [Eq. (36)].

15 else

16 Reset Bk+1←(yTs/yTy)I.

17 k←k+1.

18 until(1/2)∇f(σtk+1)TBk+1∇f(σtk+1)≤ϵ;

19 B←Bk+1, σt←σtk.

Warm start: For each inner loop, the Quasi-Newton method needs an initial guess of the inverse Hessian $B_{0}$ . Instead of fixing $B_{0}$ for every inner loop, we reuse the $B_{k}$ of the last inner loop to accelerate the convergence (shown in lines 4 and 19 in Algorithm 3).

Checking feasibility: Since the Quasi-Newton method and line search estimate without constraints, the next estimate $σ_{t}^{k + 1}$ may go beyond the constraints; in our case, each element $σ_{t}^{k + 1} [b]$ in $σ_{t}^{k + 1}$ must be inside $[0, u]$ after the step size has been determined. Therefore, in line 8, we check the feasibility of the estimate $σ_{t}^{k + 1}$ for the current step size $α_{k}$ . If it exceeds the boundary of the feasible region, we pull the estimate back into the feasible region by halving the step size. If it is still outside the feasible region, then the step size is halved again. Why do we not just set the step size so that $σ_{t}^{k + 1}$ is exactly on the boundary? The reason is the log-barrier: if $σ_{t}^{k + 1}$ is on the boundary, in other words, $σ_{t}^{k + 1} [b]$ is either 0 or $u$ , then $\log (σ_{t} [b])$ or $\log (u - σ_{t} [b])$ becomes infinite, which results in numerical instability. Therefore, the procedure described above is needed.

Checking for positive definiteness: The BFGS update rules guarantee $B_{k}$ to be positive definite if $y^{T} s > 0$ and $B ≻ 0$ are satisfied. While the latter is satisfied by giving an appropriate initial guess, the former depends on the updates at each iteration. If it is not satisfied, then the BFGS updates are no longer valid, and we reset the inverse Hessian $B_{k}$ to a scaled identity⁶³ at line 16.

4.3.1.

Jacobian

Here, we represent the Jacobian of the objective function $f_{0}$ in Eq. (29). Note that the objective function $f_{0}$ in Eq. (31) can be derived in the same manner.

We first rewrite the objective function $f_{0}$ as follows:

(43)

f_{0} = \sum_{i = 1}^{N} \sum_{j = 1}^{N} {| I_{i j} - I_{0} \sum_{k = 1}^{N_{i j}} H_{i j k} e^{- σ_{t}^{T} D_{i j k}} |}^{2}

(44)

= \sum_{i = 1}^{N} \sum_{j = 1}^{N} (I_{i j}^{2} - 2 I_{i j} I_{0} \sum_{k = 1}^{N_{i j}} H_{i j k} e^{- σ_{t}^{T} D_{i j k}} + I_{0}^{2} \sum_{k = 1}^{N_{i j}} \sum_{l = 1}^{N_{i j}} H_{i j k} e^{- σ_{t}^{T} D_{i j k}} H_{i j l} e^{- σ_{t}^{T} D_{i j l}})

(45)

= \sum_{i = 1}^{N} \sum_{j = 1}^{N} [I_{i j}^{2} - 2 I_{i j} I_{0} \sum_{k = 1}^{N_{i j}} H_{i j k} e^{- σ_{t}^{T} D_{i j k}} + I_{0}^{2} \sum_{k = 1}^{N_{i j}} \sum_{l = 1}^{N_{i j}} H_{i j k} H_{i j l} e^{- σ_{t}^{T} (D_{i j k} + D_{i j l})}],

and the gradient of

f_{0}

is given by

(46)

\frac{\partial f_{0}}{\partial σ_{t}} = \sum_{i = 1}^{N} \sum_{j = 1}^{N} [2 I_{i j} I_{0} \sum_{k = 1}^{N_{i j}} H_{i j k} e^{- σ_{t}^{T} D_{i j k}} D_{i j k} - I_{0}^{2} \sum_{k = 1}^{N_{i j}} \sum_{l = 1}^{N_{i j}} H_{i j k} H_{i j l} e^{- σ_{t}^{T} (D_{i j k} + D_{i j l})} (D_{i j k} + D_{i j l})] .

To simplify the equation, we use the following notation:

(47)

E = [\begin{matrix} e^{- σ_{t}^{T} D_{i j 1}} \\ e^{- σ_{t}^{T} D_{i j 2}} \\ ⋮ \\ e^{- σ_{t}^{T} D_{i j N_{i j}}} \end{matrix}], H = [\begin{matrix} H_{i j 1} \\ H_{i j 2} \\ ⋮ \\ H_{i j N_{i j}} \end{matrix}],

(48)

D_{i j} = [\begin{matrix} D_{i j 1} \\ D_{i j 2} \\ ⋮ \\ D_{i j N_{i j}} \end{matrix}], {\tilde{D}}_{i j} = [\begin{matrix} D_{i j 1} + D_{i j 1} & D_{i j 1} + D_{i j 2} & \dots & D_{i j 1} + D_{i j N_{i j}} \\ D_{i j 2} + D_{i j 1} & D_{i j 2} + D_{i j 2} & \dots & D_{i j 2} + D_{i j N_{i j}} \\ ⋮ & ⋮ & \dots & ⋮ \\ D_{i j N_{i j}} + D_{i j 1} & D_{i j N_{i j}} + D_{i j 2} & \dots & D_{i j N_{i j}} + D_{i j N_{i j}} \end{matrix}] .

Now, $f_{0}$ and the gradient can be represented as

(49)

f_{0} = \sum_{i = 1}^{N} \sum_{j = 1}^{N} [I_{i j}^{2} - 2 I_{i j} I_{0} E^{T} H + I_{0}^{2} {(E^{T} H)}^{2}],

(50)

\frac{\partial f_{0}}{\partial σ_{t}} = \sum_{i = 1}^{N} \sum_{j = 1}^{N} (2 I_{i j} I_{0} sum [(E \times H) \otimes D_{i j}] - I_{0}^{2} sum {[(E \times H) {(E \times H)}^{T}] \otimes {\tilde{D}}_{i j}}),

where

sum []

stands for the sum over the elements of the container [Eq. (49)] of vectors,

\times

is the element-wise product, and

\otimes

denotes the tensor product, defined as

(51)

A = [\begin{matrix} a_{11} & a_{12} & \dots & a_{1 m} \\ a_{21} & a_{22} & \dots & a_{2 m} \\ ⋮ & ⋮ & \dots & ⋮ \\ a_{n 1} & a_{n 2} & \dots & a_{n m} \end{matrix}], B = [\begin{matrix} b_{11} & b_{12} & \dots & b_{1 m} \\ b_{21} & b_{22} & \dots & b_{2 m} \\ ⋮ & ⋮ & \dots & ⋮ \\ b_{n 1} & b_{n 2} & \dots & b_{n m} \end{matrix}],

(52)

A \otimes B = [\begin{matrix} a_{11} b_{11} & a_{12} b_{12} & \dots & a_{1 m} b_{1 m} \\ a_{21} b_{21} & a_{22} b_{22} & \dots & a_{2 m} b_{2 m} \\ ⋮ & ⋮ & \dots & ⋮ \\ a_{n 1} b_{n 1} & a_{n 2} b_{n 2} & \dots & a_{n m} b_{n m} \end{matrix}] .

5. Numerical Simulations

In this section, we report the results obtained by numerical simulations using the proposed model.

The following parameters have been used in Algorithm 3: $t_{init} = 1.0$ , $μ = 1.5$ , $ϵ = 10^{- 2}$ . For the line search, the range for the step size was $α_{k} \in [0, 100]$ . For the initial guess, we used $B = I$ , $σ_{t}^{0} = 0$ . For the 2-D layered medium, the grid size was set to $N = M = 20$ with square voxels of size 1 (mm), i.e., the medium is $20 (mm) \times 20 (mm)$ , and $d A = 1$ (mm). The values of the extinction coefficients are set between 1.05 and 1.55 ( ${mm}^{- 1}$ ), and the upper bound in Eq. (30) is set to $u = 2.0$ ( ${mm}^{- 1}$ ). The parameter of the Gaussian phase function is 0.2 or 0.4, and the scattering coefficient is set to $σ_{s} = 1$ ( ${mm}^{- 1}$ ). The threshold for excluding low contribution paths is $th = 0.001$ .

The ground truth and the estimated extinction coefficients are shown in Fig. 7. The matrix plots in the top row of the figure represent five different media [from (a) to (e)] used for the simulation. Each voxel $b$ is shaded in gray according to the values of the extinction coefficient $σ_{t} [b]$ , and darker gray represents larger values of $σ_{t} [b]$ . Also, the values of $σ_{t} [b]$ are displayed at each voxel. In the same manner, the middle and bottom rows show the estimated results when the following values of the parameter of the Gaussian phase function were used: $σ^{2} = 0.2$ and 0.4. Figure 8 shows the observations $I_{i j}$ in a matrix form, from which the extinction coefficients are estimated. Each element in these plots is now an observation $I_{i j}$ . We can see observations with higher values (shown in darker shades of gray in the plots) on the diagonal. The observations obtained for $σ^{2} = 0.4$ seem to be fainter than those obtained for $σ^{2} = 0.2$ due to the larger amount of scattering.

Fig. 7

Numerical simulation results for a grid of size $20 \times 20 mm$ . Darker shades of gray represent larger values (more light is absorbed at the voxel). The bars on the side show extinction coefficient values [ ${mm}^{- 1}$ ] in gray scale. The first row shows ground truth for five different types of media [(a)–(e)] used for the simulation. The second and third rows show estimated results for $σ_{2} = 0.2$ and $σ_{2} = 0.4$ , respectively, of the Gaussian phase function.

Fig. 8

Visualization of the observations $I_{i j}$ in a matrix form. Each matrix shows $I_{i j}$ in its $i$ ’th row and $j$ ’th column. The horizontal index $i$ indicates the location of the light source, and the vertical index $j$ indicates the location of the detector. Hence, $I_{i j}$ is the light intensity with the detector at $j$ and the light source at $i$ . Darker shades of gray represent larger observation values (brighter light is observed). Left to right columns: five different media [(a)–(e)] used for the simulation in the same order as in Fig. 7. Top to bottom rows: $I_{i j}$ for T2B and L2R configurations for $σ^{2} = 0.2$ and $σ^{2} = 0.4$ .

The left-most column of Fig. 7(a) shows the simplest case: the medium has almost homogeneous extinction coefficients of value 1.05 (voxels shaded in light gray) except for a few voxels with much higher coefficients of 1.2 (voxels shaded in dark gray), which means that those voxels absorb much more light than other voxels. The coefficients are estimated reasonably well as shown in the middle and bottom rows, and the root mean squared error (RMSE) shown in Table 1 is small enough with a relative error of $0.0075 / 1.05 = 0.7 %$ to the background coefficient value. The other media, shown in columns (b) to (e), have more complex distributions of the extinction coefficients. We summarize the quality of the estimated results in terms of RMSE in Table 1. Numbers in the brackets are relative errors of RMSE to the background extinction coefficient values (i.e., 1.05). Computation time is also shown in Table 1. Note that our proposed method has been currently implemented in MATLAB®, which can be accelerated further by using C++.

Table 1

Root mean squared errors (RMSEs) and computation time for the numerical simulations for five different types of media [(a) to (e)] with a grid size of 20×20, for two different Gaussian phase function parameter values. Numbers in the brackets are relative errors of RMSE to the background extinction coefficient values (i.e., 1.05).

		(a)	(b)	(c)	(d)	(e)
RMSE	σ2=0.2	0.0067506	0.014253	0.017771	0.016220	0.057692
	σ2=0.2	(0.643%)	(1.36%)	(1.69%)	(1.54%)	(5.49%)
	σ2=0.4	0.0075305	0.014369	0.017704	0.015692	0.058464
	σ2=0.4	(0.717%)	(1.37%)	(1.69%)	(1.49%)	(5.57%)
Computation time (s)	σ2=0.2	142	113	297	190	269
Computation time (s)	σ2=0.4	127	110	186	156	267

The values of the cost function $f_{0}$ over iterations of the outer loop in Algorithm 3 are shown in Fig. 9 for each medium. These curves show that the proposed method effectively minimizes the original objective function [Eq. (31)] for the five different types of media shown here and probably for other media. Figure 10 demonstrates how the log-barriered cost function $f$ in Algorithm 3 evolves over all iterations of the inner loop; the number of iterations in the horizontal axis accumulates all inner iterations of the Quasi-Newton method. We can see that each inner loop successively minimizes the log-barriered function and the warm start (reusing the Hessian from the previous outer loop) may help the gap of values between inner loops.

Fig. 9

Original cost function values $f_{0}$ over iterations of the outer loop of Algorithm 3 with (a) $σ^{2} = 0.2$ and (b) 0.4. The horizontal axis shows the number of outer iterations, and the vertical axis represents the log of the original cost function values. Different plots indicate five different types of media [(a)–(e)] used for the simulation.

Fig. 10

Log-barriered cost function values $f$ over iterations of all inner loops of Algorithm 3 for medium (e) with $σ^{2} = 0.2$ [(a) and (b)] and 0.4 [(c) and (d)]. The horizontal axis shows the number of total inner iterations accumulated across different outer loops. The vertical axis represents the original cost function values (left) in log scale and (right) in linear scale.

5.1.

Comparison Results

We compare our method to a standard DOT with FEM (Refs. 64 and 65) using different optimization methods implemented in the Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software (EIDORS).⁶⁴^,⁶⁵ The ground truth used in this comparison is shown in the top row of Figs. 11(a)–11(e): $N = M = 24$ medium of size $24 (mm) \times 24 (mm)$ with extinction coefficient distributions almost the same as those shown in Figs. 7(a)–7(e).

Fig. 11

Numerical simulation results for a grid of size $24 \times 24 mm$ , comparing our method to diffuse optical tomography (DOT) with two solvers. Darker shades of gray represent larger values (more light is absorbed at the voxel). The bars on the side show extinction coefficient values [ ${mm}^{- 1}$ ] in gray scale. First row shows the ground truth for five different types of media [(a)–(e)] used for the simulation. Second row shows the estimated results of the proposed method. Third and fourth rows show estimated results for DOT by using Gauss-Newton (GN) and primal-dual (PD) interior point solvers.

For solving DOT by EIDORS, we used $24 \times 24 \times 24 = 1152$ triangle meshes (i.e., each voxel is divided into two triangle meshes), and for the boundary condition, we placed 16 light sources and 16 detectors at the same intervals around the medium. We chose two solvers: Gauss-Newton (GN) method and primal-dual (PD) interior point method. We used $σ_{t}^{0} = 0$ as the initial guess for both our method and EIDORS.

The results obtained by our method ( $σ^{2} = 0.4$ ) and DOT with GN and PD are shown in Fig. 11. The results obtained by the proposed method are shown in the second row, which are similar to those in the third row of Fig. 7. The third row in Fig. 11 shows the results for DOT with GN. These kind of blurred results are typical for DOT estimation due to its diffusion approximation. The last row shows results for DOT with PD, which look better than those obtained for DOT with GN, but still have a tendency of overestimating the high coefficient value areas.

We summarize RMSE values and computation time for each method in Table 2 in the same format as Table 1. RMSE values of our method are two to five times smaller than those of DOT, and this demonstrates that the proposed method can achieve much more accurate results.

The current disadvantage is its large computation cost, as our method takes up to 1000 times longer than DOT. We plan to reduce the computation cost by optimizing the code using C++ and adopting other solvers.

Table 2

RMSEs and computation time for the numerical simulations for five different types of media [(a) to (e)] with grid size of 24×24, for the proposed method and diffuse optical tomography (DOT) with two solvers. Numbers in the brackets are relative errors of RMSE to the background extinction coefficient values (i.e., 1.05).

		(a)	(b)	(c)	(d)	(e)
RMSE	Ours	0.007662	0.01244	0.026602	0.021442	0.051152
	σ2=0.4	(0.730%)	(1.18%)	(2.53%)	(2.04%)	(4.87%)
	DOT (Gauss-Newton)	0.053037	0.060597	0.7605	0.059534	0.0855
	DOT (Gauss-Newton)	(5.05%)	(5.77%)	(7.53%)	(5.67%)	(8.14%)
	DOT (primal-dual)	0.052466	0.0626	0.081081	0.066042	0.080798
	DOT (primal-dual)	(5.25%)	(5.97%)	(8.11%)	(6.60%)	(8.08%)
Computation time (s)	Ours σ2=0.4	257	217	382	306	504
Computation time (s)	DOT (Gauss-Newton)	0.397	0.390	0.407	0.404	0.453
	DOT (primal-dual)	1.11	1.09	1.14	1.08	1.15

6. Conclusion with Discussion

In this paper, we have proposed a path integral based approach to optical tomography for multiple scattering in discretized participating media. Assuming the scattering coefficients and phase function are known and uniform, the extinction coefficients at each voxel in a 2-D layered medium are estimated by using an interior point method. Numerical simulation examples are shown to demonstrate that the proposed framework works better than DOT in the simplified experimental setup, while its computation cost needs to be reduced.

There are many directions for further research, including relaxing the assumption of 2-D layered scattering model to more realistic scattering with other phase functions, using paths generated by Monte Carlo based statistical methods, extending the formulation to a full 3-D scattering model, and solving the issues mentioned below.

Limitations—stability and uniqueness: The current formulation presented in this paper estimates only the extinction coefficients; the scattering coefficients and phase function parameters are assumed to be known and uniform. This is one of the limitations of the proposed method, however, it is a common limitation of optical tomography. It is known that the scattering and absorption coefficients cannot be separated from stationary measurements of light intensity,³⁴ and the solutions are not unique. Also, given stationary measurements without angle information, the problem becomes ill-posed⁶^,⁷ and hence not stable. To overcome this limitation, we need to extend the current formulation to handle other measurements that enable stability and uniqueness, such as time-dependent, frequency-dependent, or angle-dependent measurements.

Computational cost: A large part of the computational cost of the proposed method comes from the forward model prediction [Eq. (28)], which appears in the gradient computation [Eq. (7)]. It depends on the number of paths $N_{i j}$ ; we currently use about 700 paths out of all $20^{18}$ possible paths, and for each path, we need to compute path vectors $D_{i j k}$ , $D_{i j k} + D_{i j l}$ , and factors $H_{i j k}$ . A possible acceleration is the precomputation of these variables, but this would lead to a trade-off with storage cost. Each $D_{i j k}$ has dimensions of $20 \times 20 = 400$ , each pair of $i j$ has about 700 vectors of $D_{i j k}$ , and the number of pairs $i j$ (hence observations) is $20 \times 20 = 400$ . In total, $\sim 450 MB$ memory would be required even if single precision floating numbers were used for storing all $D_{i j k}$ . Fortunately, these vectors are necessarily sparse, and we have used sparse matrices to store them. However, the increase will be linear in the number of paths $N_{i j}$ and quadratic with the grid size $\max (N, M)$ . Therefore, we plan to consider more efficient implementations.

Acknowledgments

This research is supported in part by a grant from the Japan Society for the Promotion of Science (JSPS) through the Funding Program for Next Generation World-Leading Researchers (NEXT Program) initiated by the Council for Science and Technology Policy (CSTP), and by JSPS KAKENHI Grant Number 26280061.

References

1.

S. R. Arridge and M. Schweiger, “Image reconstruction in optical tomography,” Philos. Trans. R. Soc. B 352(1354), 717–726 (1997).http://dx.doi.org/10.1098/rstb.1997.0054Google Scholar

2.

S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. 42(5), 841 (1997).http://dx.doi.org/10.1088/0031-9155/42/5/008Google Scholar

3.

J. C. Hebden, S. R. Arridge and D. T. Delpy, “Optical imaging in medicine: I. Experimental techniques,” Phys. Med. Biol. 42(5), 825 (1997).http://dx.doi.org/10.1088/0031-9155/42/5/007Google Scholar

4.

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

5.

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

6.

G. Bal, “Inverse transport theory and applications,” Inverse Probl. 25(5), 053001 (2009).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/25/5/053001Google Scholar

7.

K. Ren, “Recent developments in numerical techniques for transport-based medical imaging methods,” Commun. Comput. Phys. 8, 1–50 (2010).http://dx.doi.org/10.4208/cicp.220509.200110aGoogle Scholar

8.

A. Charette, J. Boulanger and H. K. Kim, “An overview on recent radiation transport algorithm development for optical tomography imaging,” J. Quant. Spectrosc. Radiat. Transf. 109(17–18), 2743–2766 (2008).http://dx.doi.org/10.1016/j.jqsrt.2008.06.007Google Scholar

9.

A. A. Kokhanovsky, Light Scattering Media Optics: Problems and Solutions, 3rd ed., Springer, New York (2004).Google Scholar

10.

H. W. Jensen, Realistic Image Synthesis Using Photon Mapping, AK Peters, Ltd., Natick, MA (2001).Google Scholar

11.

J. T. Kajiya, “The rendering equation,” SIGGRAPH Comput. Graph. 20, 143–150 (1986).CGRADI0097-8930http://dx.doi.org/10.1145/15886.15902Google Scholar

12.

M. Pharr and G. Humphreys, Physically Based Rendering: From Theory to Implementation, 2nd ed., Morgan Kaufmann, San Francisco, CA (2010).Google Scholar

13.

E. Cerezo et al., “A survey on participating media rendering techniques,” Vis. Comput. 21(5), 303–328 (2005).VICOE50178-2789http://dx.doi.org/10.1007/s00371-005-0287-1Google Scholar

14.

N. Kurachi, The Magic of Computer Graphics, A. K. Peters, Ltd., Natick, MA (2011).Google Scholar

15.

R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York (1965).Google Scholar

16.

J. Tessendorf, “Radiative transfer as a sum over paths,” Phys. Rev. A 35(2), 872–878 (1987).http://dx.doi.org/10.1103/PhysRevA.35.872Google Scholar

17.

M. Pauly, T. Kollig and A. Keller, “Metropolis light transport for participating media,” in Proc. of the Eurographics Workshop on Rendering Techniques, pp. 11–22, Springer-Verlag, London (2000).Google Scholar

18.

E. Veach and L. J. Guibas, “Metropolis light transport,” in Proc. of the 24th Annual Conf. on Computer Graphics and Interactive Techniques, pp. 65–76, ACM Press/Addison-Wesley Publishing Co., New York (1997).Google Scholar

19.

E. P. Lafortune and Y. D. Willems, “Rendering participating media with bidirectional path tracing,” in Proc. of the Eurographics Workshop on Rendering Techniques, pp. 91–100, Springer-Verlag, London (1996).Google Scholar

20.

E. Veach, “Robust Monte Carlo methods for light transport simulation,” PhD Thesis, Stanford University (1997).Google Scholar

21.

W. Jarosz, “Efficient Monte Carlo methods for light transport in scattering media,” PhD Thesis, University of California, San Diego (2008).Google Scholar

22.

B. Segovia, “Interactive light transport with virtual point lights,” Thèse de doctorat en informatique, Université Lyon 1 (2007).Google Scholar

23.

A. W. Kristensen, “Efficient unbiased rendering using enlightened local path sampling,” PhD Thesis, Technical University of Denmark (2011).Google Scholar

24.

M. L. Adams and E. W. Larsen, “Fast iterative methods for discrete-ordinates particle transport calculations,” Prog. Nucl. Energy 40(1), 3–159 (2002).http://dx.doi.org/10.1016/S0149-1970(01)00023-3Google Scholar

25.

K. Ren, G. Bal and A. H. Hielscher, “Frequency domain optical tomography based on the equation of radiative transfer,” SIAM J. Sci. Comput. 28, 1463–1489 (2006).SJOCE31064-8275http://dx.doi.org/10.1137/040619193Google Scholar

26.

S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, Cambridge (2004).Google Scholar

27.

T. Tamaki et al., “Multiple-scattering optical tomography with layered material,” in Int. Conf. on Signal-Image Technology Internet-Based Systems, pp. 93–99 (2013).Google Scholar

28.

B. Yuan et al., “Layered optical tomography of multiple scattering media with combined constraint optimization,” in 21st Japan-Korea Joint Workshop on Frontiers of Computer Vision, pp. 1–6 (2015).Google Scholar

29.

N. J. McCormíck, “Recent developments in inverse scattering transport methods,” Transp. Theory Stat. Phys. 13(1–2), 15–28 (1984).http://dx.doi.org/10.1080/00411458408211649Google Scholar

30.

L. Florescu, V. A. Markel and J. C. Schotland, “Single-scattering optical tomography: simultaneous reconstruction of scattering and absorption,” Phys. Rev. E 81, 016602 (2010).http://dx.doi.org/10.1103/PhysRevE.81.016602Google Scholar

31.

L. Florescu, J. C. Schotland and V. A. Markel, “Single-scattering optical tomography,” Phys. Rev. E 79, 036607 (2009).http://dx.doi.org/10.1103/PhysRevE.79.036607Google Scholar

32.

M. Schweiger, A. Gibson and S. R. Arridge, “Computational aspects of diffuse optical tomography,” Comput. Sci. Eng. 5(6), 33–41 (2003).http://dx.doi.org/10.1109/MCISE.2003.1238702Google Scholar

33.

D. A. Boas et al., “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 57–75 (2001).http://dx.doi.org/10.1109/79.962278Google Scholar

34.

A. P. Gibson, J. C. Hebden and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/50/4/R01Google Scholar

35.

L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941).http://dx.doi.org/10.1086/144246Google Scholar

36.

P. Blasi, B. Saec and C. Schlick, “A rendering algorithm for discrete volume density objects,” Comput. Graph. Forum 12(3), 201–210 (1993).Google Scholar

37.

W. M. Cornette and J. G. Shanks, “Physically reasonable analytic expression for the single-scattering phase function,” Appl. Opt. 31, 3152–3160 (1992).APOPAI0003-6935http://dx.doi.org/10.1364/AO.31.003152Google Scholar

38.

T. Nishita et al., “Display of the earth taking into account atmospheric scattering,” in Proc. of the 20th Annual Conf. on Computer Graphics and Interactive Techniques, pp. 175–182, ACM, New York (1993).Google Scholar

39.

M. Keijzer et al., “Light distributions in artery tissue: Monte Carlo simulations for finite-diameter laser beams,” Lasers Surg. Med. 9(2), 148–154 (1989).http://dx.doi.org/10.1002/lsm.1900090210Google Scholar

40.

P. Edström, “A two-phase parameter estimation method for radiative transfer problems in paper industry applications,” Inverse Probl. Sci. Eng. 16(7), 927–951 (2008).INPEEY0266-5611http://dx.doi.org/10.1080/17415970802080066Google Scholar

41.

P. Edström, “Simulation and modeling of light scattering in paper and print applications,” in Light Scattering Reviews 5, and A. A. Kokhanovsky, Ed., pp. 451–475, Springer Praxis Books, Springer, Berlin Heidelberg (2010).Google Scholar

42.

Y. Dobashi et al., “An inverse problem approach for automatically adjusting the parameters for rendering clouds using photographs,” ACM Trans. Graph. 31(6), 145 (2012).http://dx.doi.org/10.1145/2366145.2366164Google Scholar

43.

N. Joshi, C. Donner and H. W. Jensen, “Noninvasive measurement of scattering anisotropy in turbid materials by nonnormal incident illumination,” Opt. Lett. 31, 936–938 (2006).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.31.000936Google Scholar

44.

M. Baba et al., “Estimation of scattering properties of participating media using multiple-scattering renderer,” presented at Proc. of the Fourth IIEEJ Int. Workshop on Image Electronics and Visual Computing, IEEEJ, Koh Samui, Thailand (7–10 October 2014).Google Scholar

45.

I. Gkioulekas et al., “Inverse volume rendering with material dictionaries,” ACM Trans. Graph. 32, 162 (2013).ATGRDF0730-0301http://dx.doi.org/10.1145/2508363.2508377Google Scholar

46.

V. S. Antyufeev, Monte Carlo Method for Solving Inverse Problems of Radiation Transfer, Inverse and Ill-Posed Problems Series, VSP Intl. Science, The Netherlands (2000).Google Scholar

47.

G. I. Marchuk et al., “Monte Carlo methods for solving direct and inverse problems of the theory of radiative transfer in a spherical atmosphere,” in The Monte Carlo Methods in Atmospheric Optics, G. I. Marchuk et al., Eds., pp. 54–146, Springer, Berlin Heidelberg (1980).Google Scholar

48.

P. Edström, “A fast and stable solution method for the radiative transfer problem,” SIAM Rev. 47(3), 447–468 (2005).SIREAD0036-1445http://dx.doi.org/10.1137/S0036144503438718Google Scholar

49.

C. K. Hayakawa, J. Spanier and V. Venugopalan, “Coupled forward-adjoint Monte Carlo simulations of radiative transport for the study of optical probe design in heterogeneous tissues,” SIAM J. Appl. Math. 68(1), 253–270 (2007).http://dx.doi.org/10.1137/060653111Google Scholar

50.

C. Hayakawa, J. Spanier, “Perturbation Monte Carlo methods for the solution of inverse problems,” in Monte Carlo and Quasi-Monte Carlo Methods 2002, and H. Niederreiter, Ed., pp. 227–241, Springer, Berlin Heidelberg (2004).Google Scholar

51.

Y. Ishii et al., “Scattering tomography by Monte Carlo voting,” in Proc. IAPR International Conference on Machine Vision Applications (MVA2013), The International Association for Pattern Recognition (IAPR), 2013, http://www.iapr.org/ (16 July 2015).Google Scholar

52.

R. L. Barbour et al., “Model for 3-D optical imaging of tissue,” in 10th Annual Int. Geoscience and Remote Sensing Symp., Vol. 2, pp. 1395–1399 (1990).http://dx.doi.org/10.1109/IGARSS.1990.688761Google Scholar

53.

M. Raab, D. Seibert and A. Keller, “Unbiased global illumination with participating media,” in Monte Carlo and Quasi-Monte Carlo Methods, pp. 591–605 (2008).Google Scholar

54.

S. Premože, M. Ashikhmin and P. Shirley, “Path integration for light transport in volumes,” in Proc. of the 14th Eurographics Workshop on Rendering, pp. 52–63, Eurographics Association (2003).Google Scholar

55.

J. Křivánek et al., “Unifying points, beams, and paths in volumetric light transport simulation,” ACM Trans. Graph. 33(103), 1–13 (2014).ATGRDF0730-0301http://dx.doi.org/10.1145/2601097.2601219Google Scholar

56.

S. L. Jacques and X. Wang, “Path integral description of light transport in tissues,” Proc. SPIE 2979, 488–499 (1997).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.280283Google Scholar

57.

M. J. Wilson and R. K. Wang, “A path-integral model of light scattered by turbid media,” J. Phys. B 34(8), 1453 (2001).http://dx.doi.org/10.1088/0953-4075/34/8/310Google Scholar

58.

J. C. Schotland, “Tomography with diffusing photons: a Feynman path integral perspective,” in Proc. of IEEE Int. Conf. on Acoustics, Speech, and Signal Processing, Vol. 6, pp. 3791–3794 (2000).Google Scholar

59.

J. C. Schotland, “Path integrals and optical tomography,” Contemp. Math. 548, 77–84 (2011).Google Scholar

60.

P. Dutre et al., Advanced Global Illumination, 2nd ed., A K Peters/CRC Press, Boca Raton, FL (2003).Google Scholar

61.

V. S. Antyufeev, “On the distribution of a random variable,” Sib. Zh. Ind. Mat. 15(2), 3–10 (2012).Google Scholar

62.

J. Heino et al., “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).http://dx.doi.org/10.1103/PhysRevE.68.031908Google Scholar

63.

J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed., Springer Series in Operations Research and Financial Engineering, Springer, New York (2006).Google Scholar

64.

N. Polydorides, “Image reconstruction algorithms for soft-field tomography,” PhD Thesis, University of Manchester: UMIST (2002).Google Scholar

65.

N. Polydorides and W. R. Lionheart, “A MATLAB toolkit for three-dimensional electrical impedance tomography: a contribution to the electrical impedance and diffuse optical reconstruction software project,” Meas. Sci. Technol. 13(12), 1871 (2002).http://dx.doi.org/10.1088/0957-0233/13/12/310Google Scholar

Biography

Bingzhi Yuan received his BE degree in software engineering from the Beijing University of Posts and Telecommunications, China, and his ME degree in engineering from Hiroshima University, Japan, in 2010 and 2013, respectively. Currently, he is a PhD student at Hiroshima University.

Toru Tamaki received his BE, ME, and PhD degrees in information engineering from Nagoya University, Japan, in 1996, 1998, and 2001, respectively. After being an assistant professor at Niigata University, Japan, from 2001 to 2005, he is currently an associate professor in the Department of Information Engineering, Graduate School of Engineering, Hiroshima University, Japan. His research interests include computer vision and image recognition.

Takahiro Kushida received his BE degree in engineering from Hiroshima University, Japan, in 2015.

Yasuhiro Mukaigawa received his ME and PhD degrees from the University of Tsukuba in 1994 and 1997, respectively. He became a research associate at Okayama University in 1997, an assistant professor at the University of Tsukuba in 2003, an associate professor at Osaka University in 2004, and a professor at Nara Institute of Science and Technology (NAIST) in 2014. His current research interests include photometric analysis and computational photography.

Hiroyuki Kubo received his ME and PhD degrees from Waseda University in 2008 and 2012. Since 2014, he has been an assistant professor at NAIST, Japan.

Bisser Raytchev received his PhD in informatics from Tsukuba University, Japan, in 2000. After being a research associate at NTT Communication Science Labs and AIST, he is presently an assistant professor in the Department of Information Engineering, Hiroshima University, Japan. His current research interests include computer vision, pattern recognition, high-dimensional data analysis, and image processing.

Kazufumi Kaneda is a professor in the Department of Information Engineering at Hiroshima University. He received his BE, ME, and DE degrees from Hiroshima University, Japan, in 1982, 1984, and 1991, respectively. In 1986, he joined Hiroshima University. He was a visiting researcher at Brigham Young University from 1991 to 1992. His research interests include computer graphics and scientific visualization.

Citation Download Citation

Bingzhi Yuan, Toru Tamaki, Takahiro Kushida, Yasuhiro Mukaigawa, Hiroyuki Kubo, Bisser Raytchev, Kazufumi Kaneda, "Optical tomography with discretized path integral," Journal of Medical Imaging 2(3), 033501 (13 August 2015). http://dx.doi.org/10.1117/1.JMI.2.3.033501

Submission: Received ; Accepted

Access the abstract

JOURNAL ARTICLE
16 PAGES

DOWNLOAD PDF + Save for Later

SHARE

GET CITATION

< Previous Article

|

Next Article >

KEYWORDS

Scattering

Optical tomography

Light sources

Sensors

Light scattering

Light

Mass attenuation coefficient

Show All Keywords

Subscribe to Digital Library

Receive Journal Email Alerts

Create, Organize, Share

1.

Introduction

(1)

2.

Related Work

3.

Method: Forward Problem

(2)

3.1.

Forward Model

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

3.2.

Assumptions on the Path Integral Formulation

(13)

(14)

(15)

(16)

(17)

3.3.

Discretization of the Forward Model

Fig. 1

(18)

(19)

(20)

(21)

(22)

3.4.

Two-Dimensional Layered Model of Forward Scattering

Fig. 2

3.5.

Approximating the Phase Function with a Gaussian

(23)

Fig. 3

(24)

Fig. 4

Fig. 5

(25)

(26)

(27)

3.6.

Observation Model

(28)

Algorithm 1

4.

Method: Inverse Problem

4.1.

Cost Function

(29)

(30)

Fig. 6

(31)

4.2.

Optimization Problem with Inequality Constraints

4.2.1.

Unconstrained problem: Quasi-Newton

(32)

(33)

(34)

(35)

(36)

Algorithm 2

4.2.2.

Constrained problem: interior point

(37)

(38)

(39)

(40)

(41)

(42)

4.3.