19 July 2017 Primal-dual approach to optical tomography with discretized path integral with efficient formulations
Abstract
We propose an efficient optical tomography with discretized path integral. We first introduce the primal-dual approach to solve the inverse problem formulated as a constraint optimization problem. Next, we develop efficient formulations for computing Jacobian and Hessian of the cost function of the constraint nonlinear optimization problem. Numerical experiments show that the proposed formulation is faster (23.14±1.32  s) than the previous work with the log-barrier interior point method (91.17±1.48  s) for the Shepp–Logan phantom with a grid size of 24×24, while keeping the quality of the estimation results (root-mean-square error increasing by up to 12%).

## Introduction

Optical tomography18 is a technique for reconstructing inside an object by illuminating it with a light probe and observing the light penetrating through the object. In contrast to x-ray computed tomography, which uses x-rays instead of light, safer tomographic methods are demanded and scattering optical tomography methods are recently attracting computer vision researchers’ attention.915

We investigate a recently proposed optical tomography with discretized path integral.12,13,15 Their method takes advantages of the path integral formulation and formulates the inverse problem of optical tomography as a constraint nonlinear least squared optimization problem. This method benefits from various optimization techniques, which is not the case for voting11,14 or genetic algorithms.10 They solved the constraint optimization problem by using the log-barrier (LB) interior point method16 with inner loops of Newton method12 and quasi-Newton method.15 They have shown that optical tomography with discretized path integral produces better estimation results compared to a standard diffuse optical tomography (DOT),15 however, its high computation cost is a problem for further development.

In this paper, we propose two contributions to tackle the problem. First, we introduce the primal-dual (PD) approach to solve the constraint optimization because it is known that the PD interior point method is more efficient than the LB method.16 Second, we propose new formulations of equations. The main bottle-neck of the previous approaches12,13,15 is Jacobian and Hessian which are computationally demanding as the number of paths increases. Our new formulations are equivalent to the previous ones, but much more efficient. Numerical simulations show that the proposed approach accelerates the estimation by a factor of 100. (Conference versions of this paper were presented.17,18 This paper extends those versions with the extended description of the PD approach and the efficient formulations, and additional experiments with optimized codes with comparisons.)

There exists a number of acceleration methods for tomographic reconstruction; such as dimension reduction,19 power acceleration,20 and, most importantly, graphic processing units (GPUs).21 GPUs have been used for accelerating a forward problem22 and forward and backward projections.2324.25 In addition, thanks to the recent progress of general-purpose computing on GPUs, GPUs have become popular for speeding up iterative computations for solving compressive sensing formulations26 and large linear systems.27,28 We have not implemented the proposed approach with any GPUs in this paper, however, the use of GPUs would be beneficial for the proposed formulations because computing a Jacobian matrix can be further accelerated.29

In the following, we briefly summarize the constraint problem to be solved in Sec. 2. Then in Sec. 3, we develop the PD method to solve the problem by taking into account the structure of the problem. In Sec. 4, we derive new efficient formulations to compute Jacobian and Hessian, with comparisons to and discussions of the previous formulations. We show simulation results in Sec. 5 to show the improvement of the proposed method in terms of computation cost.

## Preliminary

In this section, we briefly review the formulation of optical tomography with a discretized path integral. Details can be found in Ref. 15.

We are interested in estimating the extinction coefficients ${\mathbf{\sigma }}_{t}$ of each voxel in a two-dimensional (2-D) medium as shown in Fig. 1. Extinction coefficients represent how much light attenuates at each point. We follow the 2-D layer model,15 that is, we assume the following layer scattering with the following properties. Suppose that the 2-D medium has a layered structure and is discretized into voxels of an $M×N$ grid; it has $M$ layers of $N$ voxels. Therefore, the problem is to estimate extinction coefficients of each voxel in the grid.

## Fig. 1

2-D layered model of scattering.15 This example of a $5×5$ grid shows a path 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.

With this layered medium, we use 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 (see Fig. 2). More specifically, the light source point ${x}_{0}$ is located on the boundary of the top surface of the voxels in the top layer. The detector point ${x}_{M+1}$ is located on the boundary of the bottom surface of the voxels in the bottom layer. Then, forward scattering happens layer by layer; light is scattered at the center of a voxel in a layer and goes to the center of a voxel in the next (below) layer. By connecting the centers ${x}_{1},\dots ,{x}_{M}$ of voxels of each layer, we have a path ${x}_{0},{x}_{1},\dots ,{x}_{M},{x}_{M+1}$ of light scattering connecting the light source and the detector. Let $i$ and $j$ be voxel indices of the light source and detector locations, respectively. By restricting the light paths only to those connecting $i$ and $j$, the observed light ${I}_{ij}$ by the detector is the sum of contributions of all paths connecting $i$ and $j$. This is written as follows:

## (1)

${I}_{ij}={I}_{0}\sum _{k=1}^{{N}_{ij}}{H}_{ijk}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}},$
where ${\mathbf{\sigma }}_{t}\in {\mathbb{R}}^{NM}$ is a vector to be estimated, and each element is the extinction coefficient of a voxel, as shown Fig. 2. Vector ${\mathbf{D}}_{k}\in {\mathbb{R}}^{NM}$ represents a complete path $k$ connecting $i$ and $j$, and each element is the length of the part of the path segment passing through the corresponding voxel. Factor ${H}_{ijk}$ encodes scattering coefficients and the phase function, ${I}_{0}$ is the intensity of the light source, and ${N}_{ij}$ is the number of paths connecting $i$ and $j$. Parameters ${H}_{ijk}$, ${I}_{0}$, ${N}_{ij}$, ${\mathbf{D}}_{k}$ are given and fixed. The problem is to estimate ${\mathbf{\sigma }}_{t}$ based on observations ${I}_{ij}$.

## Fig. 2

An example of a path in a $5×5$ grid. The path ${\mathbf{D}}_{ijk}$ is one of paths that connect locations $i$ and $j$. The grid is serialized to vector ${\mathbf{\sigma }}_{t}$, and also to vector ${\mathbf{D}}_{ijk}$ separately, but in a consistent order. These vectors are used to represent the exponential attenuation of light along the path by inner product followed by the exponential function as $\mathrm{exp}\left(-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}\right)$ in Eq. (3).

By changing positions of the light source $i$ and the detector $j$, we obtain a set of observations $\left\{{I}_{ij}\right\}$, resulting in the following nonlinear least squares problem under box constraints for extinction coefficients to be positive

## (2)

$\underset{{\mathbf{\sigma }}_{t}}{\mathrm{min}}\text{\hspace{0.17em}\hspace{0.17em}}f\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{s.t.}\text{\hspace{0.17em}\hspace{0.17em}}{\sigma }_{tl}\le {\mathbf{\sigma }}_{t}\le {\sigma }_{tu},$
where $f$ is the cost function

## (3)

$f=\sum _{i=1}^{N}\sum _{j=1}^{N}{|{I}_{ij}-{I}_{0}\sum _{k=1}^{{N}_{ij}}{H}_{ijk}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}}|}^{2},$
and ${\sigma }_{tl}$ and ${\sigma }_{tu}$ are lower and upper bounds, respectively. The box constraints are due to the nature of the extinction coefficient being positive (i.e., ${\sigma }_{tl}>0$) and the numerical stability of excluding unrealistic large values.

## Primal-Dual Interior Point Method of the Inverse Problem

Here, we develop a PD interior point method to solve the inverse problem [Eq. (2)]. It is an inequality constraint optimization with box constraints, which is straightforward in applying a standard PD method.16 However, we can use the structure of the box constraints, hence, we derive an efficient algorithm by using the problem structure.

## 3.1.

### Primal-Dual Method

We first rewrite the inequality constraint problem to an equivalent problem with equality constraints with slack variables $\mathbf{s}={\left({s}_{1},\dots ,{s}_{2NM}\right)}^{\mathrm{T}}$ as follows:

## (4)

$\underset{{\mathbf{\sigma }}_{t}}{\mathrm{min}}\text{\hspace{0.17em}\hspace{0.17em}}f\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{s.t.}\text{\hspace{0.17em}\hspace{0.17em}}\mathbf{c}-\mathbf{s}=\mathbf{0},\phantom{\rule[-0.0ex]{1em}{0.0ex}}0\le \mathbf{s},$
where $\mathbf{c}$ is a vector of the box constraints

## (5)

$\mathbf{c}=\left(\begin{array}{c}{c}_{1}\left({\mathbf{\sigma }}_{t}\right)\\ ⋮\\ {c}_{2NM}\left({\mathbf{\sigma }}_{t}\right)\end{array}\right)=\left(\begin{array}{c}{\mathbf{\sigma }}_{t}-{\sigma }_{tl}\mathbf{1}\\ {\sigma }_{tu}\mathbf{1}-{\mathbf{\sigma }}_{t}\end{array}\right).$
Here, ${c}_{i}$ is the $i$’th constraint, and $\mathbf{1}$ is a vector of ones.

The Lagrangian $L$ of the above problem is

## (6)

$L\left({\mathbf{\sigma }}_{t},\mathbf{s},\mathbf{z}\right)=f-\sum _{i=1}^{2NM}{z}_{i}\left({c}_{i}-{s}_{i}\right)=f-{\mathbf{z}}^{\mathrm{T}}\left(\mathbf{c}-\mathbf{s}\right),$
where $\mathbf{z}$ is a vector of Lagrangian multipliers or dual variables.

The KKT conditions of Eq. (4) with duality gap $\mu$ is written as

## (7)

$\nabla f-{A}^{\mathrm{T}}\mathbf{z}=\mathbf{0},\phantom{\rule[-0.0ex]{1em}{0.0ex}}S\mathbf{z}-\mu \mathbf{1}=\mathbf{0},\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathbf{c}-\mathbf{s}=\mathbf{0},\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathbf{s}\ge \mathbf{0},\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathbf{z}\ge \mathbf{0},$
where $S=\mathrm{diag}\left(\mathbf{s}\right)$, and

## (8)

$A=\nabla \mathbf{c}=\left(\begin{array}{c}\nabla {c}_{1}\\ ⋮\\ \nabla {c}_{2NM}\end{array}\right)=\left(\begin{array}{c}I\\ -I\end{array}\right).$
Here, $I$ is an identity matrix.

To solve the system of the KKT conditions by using Newton’s method, we have the following system of equations:

## (9)

$\left(\begin{array}{ccc}{\nabla }^{2}L& \mathbf{0}& -{A}^{\mathrm{T}}\\ \mathbf{0}& Z& S\\ A& -I& \mathbf{0}\end{array}\right)\left(\begin{array}{c}{\mathbf{p}}_{\sigma }\\ {\mathbf{p}}_{s}\\ {\mathbf{p}}_{z}\end{array}\right)=-\left(\begin{array}{c}\nabla f-{A}^{\mathrm{T}}\mathbf{z}\\ S\mathbf{z}-\mu \mathbf{1}\\ \mathbf{c}-\mathbf{s}\end{array}\right),$
where $Z=\mathrm{diag}\left(\mathbf{z}\right)$ and ${\nabla }^{2}L$ is the Hessian of the Lagrangian.

## 3.2.

### Solving the System Efficiently

The matrix in Eq. (9) is of the size $5NM×5NM$, which is sparse but large and computationally expensive to solve. We, therefore, develop an efficient way to solve the system by using the problem structure.

First, the system is explicitly written as follows:

## (10)

$\left\{\begin{array}{l}{\nabla }^{2}L{\mathbf{p}}_{\sigma }-{A}^{\mathrm{T}}{\mathbf{p}}_{z}=-\nabla f+{A}^{\mathrm{T}}\mathbf{z}\\ Z{\mathbf{p}}_{s}+S{\mathbf{p}}_{z}=-S\mathbf{z}+\mu \mathbf{1}\\ A{\mathbf{p}}_{\sigma }-{\mathbf{p}}_{s}=-\mathbf{c}+\mathbf{s}\end{array}.$
Substituting the last equation into the second one yields

## (11)

${A}^{\mathrm{T}}{S}^{-1}Z\left(A{\mathbf{p}}_{\sigma }+\mathbf{c}-\mathbf{s}\right)+{A}^{\mathrm{T}}{\mathbf{p}}_{z}=-{A}^{\mathrm{T}}\mathbf{z}+\mu {A}^{\mathrm{T}}{S}^{-1}\mathbf{1},$
and we add both sides to the first equation to obtain

## (12)

${\nabla }^{2}L{\mathbf{p}}_{\sigma }+{A}^{\mathrm{T}}{S}^{-1}Z\left(A{\mathbf{p}}_{\sigma }+\mathbf{c}-\mathbf{s}\right)=-{A}^{\mathrm{T}}\mathbf{z}+\mu {A}^{\mathrm{T}}{S}^{-1}\mathbf{1}-\nabla f+{A}^{\mathrm{T}}\mathbf{z}.$
Here, we define $\mathbf{w}={\mathbf{s}}^{-1}\odot \mathbf{z}$ and $\mathbf{y}=\mu {\mathbf{s}}^{-1}-\mathbf{w}\odot \mathbf{c}+\mathbf{z}$, where $\odot$ is the Hadamard (element-wise) product, and ${\mathbf{s}}^{-1}$ is a vector of element-wise reciprocals of $\mathbf{s}$. Then, we have

## (13)

${\mathbf{p}}_{\sigma }={\left({\nabla }^{2}L+{A}^{\mathrm{T}}{S}^{-1}ZA\right)}^{-1}\left(-\nabla f+{A}^{\mathrm{T}}\mathbf{y}\right),$

## (14)

${\mathbf{p}}_{z}=\mu {\mathbf{s}}^{-1}-\mathbf{w}\odot {\mathbf{p}}_{s}-\mathbf{z}.$
By exploiting the structure of matrix $A$ and defining $\mathbf{w}={\left({\mathbf{w}}_{l}^{\mathrm{T}},{\mathbf{w}}_{u}^{\mathrm{T}}\right)}^{\mathrm{T}}$ to split a vector into two parts corresponding to lower and upper bound constraints, we have

## (15)

${A}^{\mathrm{T}}{S}^{-1}ZA=\left(\begin{array}{cc}I& -I\end{array}\right)\mathrm{diag}\left(\mathbf{w}\right)\left(\begin{array}{c}I\\ -I\end{array}\right)=\mathrm{diag}\left({\mathbf{w}}_{l}+{\mathbf{w}}_{u}\right).$
Similarly, we define $\mathbf{y}={\left({\mathbf{y}}_{l}^{\mathrm{T}},{\mathbf{y}}_{u}^{\mathrm{T}}\right)}^{\mathrm{T}}$ to simplify ${A}^{\mathrm{T}}\mathbf{y}$ as ${A}^{\mathrm{T}}\mathbf{y}={\mathbf{y}}_{l}-{\mathbf{y}}_{u}$ and then

## (16)

$A{\mathbf{p}}_{\sigma }=\left(\begin{array}{c}I\\ -I\end{array}\right){\mathbf{p}}_{\sigma }=\left(\begin{array}{c}{\mathbf{p}}_{\sigma }\\ -{\mathbf{p}}_{\sigma }\end{array}\right).$
Now, the solution is given by

## (17)

$\left\{\begin{array}{l}{\mathbf{p}}_{\sigma }={\left[{\nabla }^{2}L+\mathrm{diag}\left({\mathbf{w}}_{l}+{\mathbf{w}}_{u}\right)\right]}^{-1}\left(-\nabla f+{\mathbf{y}}_{l}-{\mathbf{y}}_{u}\right)\\ {\mathbf{p}}_{s}=\left(\begin{array}{c}{\mathbf{p}}_{\sigma }\\ -{\mathbf{p}}_{\sigma }\end{array}\right)+\mathbf{c}-\mathbf{s}\\ {\mathbf{p}}_{z}=\mu {\mathbf{s}}^{-1}-\mathbf{w}\odot {\mathbf{p}}_{s}-\mathbf{z}\end{array},$
which involves the inversion of the size $NM×NM$.

## 3.3.

### Update Variables

Once ${\mathbf{p}}_{\sigma }$, ${\mathbf{p}}_{s}$, and ${\mathbf{p}}_{z}$ are obtained, we then estimate the step length to update the parameters.16 The maximum of the step lengths is given by the following rule:

## (18)

$\left\{\begin{array}{c}{\alpha }_{s}^{\mathrm{max}}=\mathrm{max}\left\{\alpha \in \left[0,1\right]|\mathbf{s}+\alpha {\mathbf{p}}_{s}\ge \left(1-\tau \right)\mathbf{s}\right\}\\ {\alpha }_{z}^{\mathrm{max}}=\mathrm{max}\left\{\alpha \in \left[0,1\right]|\mathbf{z}+\alpha {\mathbf{p}}_{z}\ge \left(1-\tau \right)\mathbf{z}\right\}\end{array},$
with $\tau \in \left(0,1\right)$ used (for example, $\tau =0.995$). This prevents variables $\mathbf{s}$ and $\mathbf{z}$ from approaching the lower boundary.

Next, we perform the backtracking line search30 to estimate acceptable step lengths ${\alpha }_{s}$ and ${\alpha }_{z}$. To this end, we use the following exact merit function $\varphi$ with $\eta \in \left(0,1\right)$:

## (19)

$\varphi \left({\mathbf{\sigma }}_{t},\mathbf{s}\right)=f-\mu \sum _{i=1}^{2MN}\mathrm{log}\left({s}_{i}\right)+v‖\mathbf{c}\left(x\right)-\mathbf{s}‖,$
and make a sufficient decrease requirement

## (20)

$\varphi \left({\mathbf{\sigma }}_{t}+{\alpha }_{s}{\mathbf{p}}_{\sigma },\mathbf{s}+{\alpha }_{s}{\mathbf{p}}_{s}\right)\le \varphi \left({\mathbf{\sigma }}_{t},\mathbf{s}\right)+\eta {\alpha }_{s}{D}_{\varphi }\left({\mathbf{\sigma }}_{t},\mathbf{s};{\mathbf{p}}_{\sigma },{\mathbf{p}}_{s}\right),$
where ${D}_{\varphi }\left({\mathbf{\sigma }}_{t},\mathbf{s};{\mathbf{p}}_{\sigma },{\mathbf{p}}_{s}\right)$ denotes the directional derivative of $\varphi$ in the direction $\left({\mathbf{p}}_{\sigma },{\mathbf{p}}_{s}\right)$. The step lengths ${\alpha }_{s}$ and ${\alpha }_{z}$ are found in the ranges ${\alpha }_{s}\in \left(0,{\alpha }_{s}^{\mathrm{max}}\right]$ and ${\alpha }_{z}\in \left(0,{\alpha }_{z}^{\mathrm{max}}\right]$ so that Eq. (20) is satisfied.

Then, the parameters ${\mathbf{\sigma }}_{t}$, $\mathbf{s}$, and $\mathbf{z}$ are updated as

## (21)

$\left\{\begin{array}{l}{\mathbf{\sigma }}_{t}←{\mathbf{\sigma }}_{t}+{\alpha }_{s}{\mathbf{p}}_{\sigma }\\ \mathbf{s}←\mathbf{s}+{\alpha }_{s}{\mathbf{p}}_{s}\\ \mathbf{z}←\mathbf{z}+{\alpha }_{z}{\mathbf{p}}_{z}\end{array}.$

Once the following error function is smaller than a given threshold, the PD interior point method stops

## (22)

$E\left({\mathbf{\sigma }}_{t},\mathbf{s},\mathbf{z};\mu \right)=\mathrm{max}\left\{‖\nabla f-{A}^{\mathrm{T}}\mathbf{z}‖,‖S\mathbf{z}-\mu \mathbf{1}‖,‖\mathbf{c}-\mathbf{s}‖\right\}.$

Algorithm 1 summarizes the PD interior point method developed above. Note that the Hessian ${\nabla }^{2}L$ can be approximated as ${B}_{k}$ at each iteration $k$ by the quasi-Newton method, instead of the full Hessian used by Newton’s method. We will compare Newton’s method and the quasi-Newton method in the section of experiments.

## Algorithm 1

Primal-dual interior point with line search.

 Data:$\mu >0$, $\sigma \in \left(0,1\right)$, ${\epsilon }_{\mathrm{TOL}}>0$, ${\epsilon }_{\mu }$, $\eta \in \left(0,1\right)$, $k=0$. Input: A feasible initial estimates ${\sigma }_{t}$, $\mathbf{s}>0$, and $\mathbf{z}>0$. Input:${B}_{0}=I$ // Only for quasi-Newton Result: Estimates of ${\sigma }_{t}$ 1 repeat 2 repeat // inner loop 3  Compute the decent direction $\mathbf{p}=\left({\mathbf{p}}_{\sigma },{\mathbf{p}}_{s},{\mathbf{p}}_{z}\right)$ 4  Compute the step lengths ${\alpha }_{s}$ and ${\alpha }_{z}$ 5  Update ${\mathbf{\sigma }}_{t}$, $\mathbf{s}$, $\mathbf{z}$ 6  Update the approximation ${B}_{k}$ // Only for quasi-Newton 7  Set $k←k+1$ 8 until$E\left({\mathbf{\sigma }}_{tk},{\mathbf{s}}_{k},{\mathbf{z}}_{k};\mu \right)\le {\epsilon }_{\mu }$; 9 $\mu ←\sigma \mu$ 10 ${\epsilon }_{\mu }←\mu$ 11 until$E\left({\mathbf{\sigma }}_{tk},{\mathbf{s}}_{k},{\mathbf{z}}_{k};0\right)\le {\epsilon }_{\mathrm{TOL}}$;

## Efficient Formulations

The most computationally intensive part of the PD algorithm shown above is the computation of Hessians for Newton’s method and Jacobians for Newton’s and quasi-Newton methods. We propose here efficient formulations of Hessian and Jacobian of the problem, whose computational cost is much smaller than naive formulations used in the previous approaches.

First, we show the naive and old formulations of Hessian and Jacobian, then introduce our new formulations, followed by discussions on computational cost.

## 4.1.

### Previous Old Formulations for Inverse Problem

Here, we show how the previous approaches12,13,15 do. We call these the old formulations.

## 4.1.1.

#### Jacobian: old formulation

Remember that the objective function $f$ to be minimized is

## (23)

$f=\sum _{i=1}^{N}\sum _{j=1}^{N}\left[{I}_{ij}^{2}-2{I}_{ij}{I}_{0}\sum _{k=1}^{{N}_{ij}}{H}_{ijk}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}}+{I}_{0}^{2}\sum _{k=1}^{{N}_{ij}}\sum _{l=1}^{{N}_{ij}}{H}_{ijk}{H}_{ijl}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}\left({\mathbf{D}}_{ijk}+{\mathbf{D}}_{ijl}\right)}\right].$
The gradient of $f$ is given as follows by taking the derivative of the objective function:

## (24)

$\frac{\partial f}{\partial {\mathbf{\sigma }}_{t}}=\sum _{i=1}^{N}\sum _{j=1}^{N}\left[2{I}_{ij}{I}_{0}\sum _{k=1}^{{N}_{ij}}{H}_{ijk}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}}{\mathbf{D}}_{ijk}\phantom{\rule{0ex}{0ex}}-{I}_{0}^{2}\sum _{k=1}^{{N}_{ij}}\sum _{l=1}^{{N}_{ij}}{H}_{ijk}{H}_{ijl}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}\left({\mathbf{D}}_{ijk}+{\mathbf{D}}_{ijl}\right)}\left({\mathbf{D}}_{ijk}+{\mathbf{D}}_{ijl}\right)\right].$
To simplify the equation, the following notations are introduced:

## (25)

${\mathbf{E}}_{ij}=\left(\begin{array}{c}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ij1}}\\ {e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ij2}}\\ ⋮\\ {e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ij{N}_{ij}}}\end{array}\right),\phantom{\rule[-0.0ex]{1em}{0.0ex}}{\mathbf{H}}_{ij}=\left(\begin{array}{c}{H}_{ij1}\\ {H}_{ij2}\\ ⋮\\ {H}_{ij{N}_{ij}}\end{array}\right),\phantom{\rule[-0.0ex]{1em}{0.0ex}}{D}_{ij}=\left(\begin{array}{c}{\mathbf{D}}_{ij1}\\ {\mathbf{D}}_{ij2}\\ ⋮\\ {\mathbf{D}}_{ij{N}_{ij}}\end{array}\right),$

## (26)

$\stackrel{˜}{{D}_{ij}}=\left(\begin{array}{ccc}{\mathbf{D}}_{ij1}+{\mathbf{D}}_{ij1}& \cdots & {\mathbf{D}}_{ij1}+{\mathbf{D}}_{ij{N}_{ij}}\\ {\mathbf{D}}_{ij2}+{\mathbf{D}}_{ij1}& \cdots & {\mathbf{D}}_{ij2}+{\mathbf{D}}_{ij{N}_{ij}}\\ ⋮& \ddots & ⋮\\ {\mathbf{D}}_{ij{N}_{ij}}+{\mathbf{D}}_{ij1}& \cdots & {\mathbf{D}}_{ij{N}_{ij}}+{\mathbf{D}}_{ij{N}_{ij}}\end{array}\right).$
Now, $f$ and the gradient are rewritten as follows:

## (27)

$f=\sum _{i=1}^{N}\sum _{j=1}^{N}\left[{I}_{ij}^{2}-2{I}_{ij}{I}_{0}{\mathbf{E}}_{ij}^{\mathrm{T}}{\mathbf{H}}_{ij}+{I}_{0}^{2}{\left({\mathbf{E}}_{ij}^{\mathrm{T}}{\mathbf{H}}_{ij}\right)}^{2}\right],$

## (28)

$\frac{\partial f}{\partial {\sigma }_{t}}=\sum _{i=1}^{N}\sum _{j=1}^{N}\left(2{I}_{ij}{I}_{0}\text{sum}\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\otimes {D}_{ij}\right]-{I}_{0}^{2}\text{sum}\left\{\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\right]\otimes \stackrel{˜}{{D}_{ij}}\right\}\right),$
where $\text{sum}\left\{\right\}$ stands for the sum over the elements of the container [Eq. (26)] of vectors, and $\otimes$ denotes the tensor product.

## 4.1.2.

#### Hessian: old formulation

We define the following notations:

## (29)

${\mathbf{\beta }}_{ijx}=\left(\begin{array}{c}{D}_{ij1x}\\ {D}_{ij2x}\\ ⋮\\ {D}_{ijkx}\end{array}\right),$

## (30)

${\mathbf{\gamma }}_{ijx}=\left(\begin{array}{ccc}{D}_{ij1x}+{D}_{ij1x}& \cdots & {D}_{ij1x}+{D}_{ijkx}\\ {D}_{ij2x}+{D}_{ij1x}& \cdots & {D}_{ij2x}+{D}_{ijkx}\\ ⋮& \ddots & ⋮\\ {D}_{ijkx}+{D}_{ij1x}& \cdots & {D}_{ijkx}+{D}_{ijkx}\end{array}\right),$
where ${D}_{ijkx}$ stands for the $x$’th element of ${\mathbf{D}}_{ijk}$.

Now, the second-order derivate of $f$ can be represented as follows:

## (31)

$\frac{{\partial }^{2}f}{{\partial }^{2}{\mathbf{\sigma }}_{t}}=\left(\begin{array}{cccc}{h}_{1,1}& {h}_{1,2}& \cdots & {h}_{1,NM}\\ {h}_{2,1}& {h}_{2,2}& \cdots & {h}_{2,NM}\\ ⋮& ⋮& \cdots & ⋮\\ {h}_{NM,1}& {h}_{NM,2}& \cdots & {h}_{NM,NM}\end{array}\right),$
where each element is given by

## (32)

${h}_{xy}=-\sum _{i=1}^{N}\sum _{j=1}^{N}{I}_{ij}{I}_{0}\text{sum}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\odot {\mathbf{\beta }}_{ijx}\odot {\mathbf{\beta }}_{ijy}\right)+\sum _{i=1}^{N}\sum _{j=1}^{N}{I}_{0}^{2}\text{sum}\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\odot {\mathbf{\gamma }}_{ijx}\odot {\mathbf{\gamma }}_{ijy}\right].$

## 4.2.

### Proposed New Efficient Formulation

The problem of the previous old formulations of Jacobian and Hessian shown above is the computation cost increasing as the number ${N}_{ij}$ of paths increases. As discussed later, the computation cost is $O\left({N}_{ij}^{2}\right)$ on average.

The idea of the proposed formulation is to explore the property of the exponential function and its derivative in the problem. As shown below, the computation cost can be reduced to $O\left({N}_{ij}\right)$ on average.

## 4.2.1.

#### Jacobian: new formulation

First, we rewrite the cost function as follows:

## (33)

$f=\sum _{i=1}^{N}\sum _{j=1}^{N}{r}_{ij}^{2},$
where ${r}_{ij}$ is a residual

## (34)

${r}_{ij}={I}_{ij}-{I}_{0}\sum _{k=1}^{{N}_{ij}}{H}_{ijk}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}}={I}_{ij}-{I}_{0}{\mathbf{E}}_{ij}^{\mathrm{T}}{\mathbf{H}}_{ij}.$

Now we use the chain rule of differentiation, and we have

## (35)

$\frac{\partial f}{\partial {\mathbf{\sigma }}_{t}}=\sum _{i=1}^{N}\sum _{j=1}^{N}2{r}_{ij}\frac{\partial {r}_{ij}}{\partial {\mathbf{\sigma }}_{t}},$
where

## (36)

$\frac{\partial {r}_{ij}}{\partial {\mathbf{\sigma }}_{t}}={I}_{0}\sum _{k=1}^{{N}_{ij}}{H}_{ijk}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}}{\mathbf{D}}_{ijk}={I}_{0}{D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right).$

Here, we define [Note that this is not the same as the one defined in the previous approaches above, which is a structure used in MATLAB codes to store ${\mathbf{D}}_{ijk}$. Here, ${D}_{ij}$ is an $\left(NM\right)×{N}_{ij}$ matrix.]

## (37)

${D}_{ij}=\left(\begin{array}{c}{\mathbf{D}}_{ij1},{\mathbf{D}}_{ij2},\dots ,{\mathbf{D}}_{ij{N}_{ij}}\end{array}\right),$
which has ${N}_{ij}$ vectors of dimension $NM$, and $\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)$ is a coefficient vector. Therefore,

## (38)

$\frac{\partial f}{\partial {\mathbf{\sigma }}_{t}}=\sum _{i=1}^{N}\sum _{j=1}^{N}2{I}_{0}\left({I}_{ij}-{I}_{0}{\mathbf{E}}_{ij}^{\mathrm{T}}{\mathbf{H}}_{ij}\right){D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right).$

## 4.2.2.

#### Discussion

Suppose that the expectation of the number of paths is $\overline{N}=E\left[{N}_{ij}\right]$. Then, the computations for the proposed new formulation of Jacobian above are:

• $O\left(\overline{N}\right)$ multiplications for ${\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}$,

• $O\left(\overline{N}\right)$ additions for ${\mathbf{E}}_{ij}^{\mathrm{T}}{\mathbf{H}}_{ij}$,

• $O\left(\overline{N}NM\right)$ multiplications and additions for ${D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)$ because there are $O\left(\overline{N}\right)$ vectors of dimension $NM$,

for each $i$ and $j$. In total, it takes $O\left(\overline{N}{N}^{3}M\right)$ operations to compute $NM$ elements of the Jacobian or $O\left(\overline{N}{N}^{2}\right)$ operations per element.

Contrary, for each $i$ and $j$, the previous old formulation Eq. (28) needs:

• $O\left(\stackrel{‾}{N}\right)$ multiplications for ${\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}$,

• $O\left(\stackrel{‾}{N}NM\right)$ multiplications for $\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\otimes {D}_{ij}$ because there are $O\left(\overline{N}\right)$ vectors of dimension $NM$,

• $O\left(\overline{N}NM\right)$ additions for $\text{sum}\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\otimes {D}_{ij}\right]$,

for the first term, and

• $O\left({\stackrel{‾}{N}}^{2}\right)$ multiplications for $\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}$,

• $O\left({\overline{N}}^{2}NM\right)$ additions for computing $\stackrel{˜}{{D}_{ij}}$ because there are $O\left({\stackrel{‾}{N}}^{2}\right)$ vectors of dimension $NM$,

• $O\left({\stackrel{‾}{N}}^{2}NM\right)$ multiplications for $\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\right]\otimes \stackrel{˜}{{D}_{ij}}$,

• $O\left({\overline{N}}^{2}NM\right)$ additions for $\text{sum}\left\{\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\right]\otimes \stackrel{˜}{{D}_{ij}}\right\}$,

for the second term. In total, it takes $O\left({\overline{N}}^{2}{N}^{3}M\right)$ operations to compute $NM$ elements of the Jacobian, or $O\left({\overline{N}}^{2}{N}^{2}\right)$ operations per element. The difference is mainly caused by the second term of Eq. (28).

In summary, the proposed new formulation has the cost of $O\left(\stackrel{‾}{N}{N}^{2}\right)$ operations per element, whereas the previous old formulation has the cost of $O\left({\overline{N}}^{2}{N}^{2}\right)$ operations per element. Table 1 summarizes the discussion above.

## Table 1

Comparison of the new and old formulations for computing the Jacobian.

TermsNewOld
${\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}$$O\left(\overline{N}\right)$$O\left(\stackrel{‾}{N}\right)$
${\mathbf{E}}_{ij}^{\mathrm{T}}{\mathbf{H}}_{ij}$$O\left(\stackrel{‾}{N}\right)$
${D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)$$O\left(\overline{N}NM\right)$
$\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\odot {D}_{ij}$$O\left(\overline{N}NM\right)$
$\text{sum}\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\odot {D}_{ij}\right]$$O\left(\stackrel{‾}{N}NM\right)$
$\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}$$O\left({\stackrel{‾}{N}}^{2}\right)$
$\stackrel{˜}{{D}_{ij}}$$O\left({\stackrel{‾}{N}}^{2}NM\right)$
$\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\right]\odot \stackrel{˜}{{D}_{ij}}$$O\left({\stackrel{‾}{N}}^{2}NM\right)$
$\text{sum}\left\{\left[\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\right]\odot \stackrel{˜}{{D}_{ij}}\right\}$$O\left({\overline{N}}^{2}NM\right)$
Per element$O\left(\stackrel{‾}{N}{N}^{2}\right)$$O\left({\overline{N}}^{2}{N}^{2}\right)$

## 4.2.3.

#### Hessian: new formulation

In the same manner, we can derive the Hessian as follows. From the Jacobian

## (39)

$\frac{\partial f}{\partial {\mathbf{\sigma }}_{t}}=\sum _{i=1}^{N}\sum _{j=1}^{N}2{r}_{ij}\frac{\partial {r}_{ij}}{\partial {\mathbf{\sigma }}_{t}},$
we obtain the Hessian by using the chain rule of differentiation

## (40)

$\frac{{\partial }^{2}f}{\partial {\mathbf{\sigma }}_{t}^{2}}=\sum _{i=1}^{N}\sum _{j=1}^{N}2\frac{\partial {r}_{ij}}{\partial {\mathbf{\sigma }}_{t}}{\frac{\partial {r}_{ij}}{\partial {\mathbf{\sigma }}_{t}}}^{\mathrm{T}}+2{r}_{ij}\frac{{\partial }^{2}{r}_{ij}}{\partial {\mathbf{\sigma }}_{t}^{2}},$
where

## (41)

$\frac{{\partial }^{2}{r}_{ij}}{\partial {\mathbf{\sigma }}_{t}^{2}}=-{I}_{0}\sum _{k=1}^{{N}_{ij}}{H}_{ijk}{e}^{-{\mathbf{\sigma }}_{t}^{\mathrm{T}}{\mathbf{D}}_{ijk}}{\mathbf{D}}_{ijk}{\mathbf{D}}_{ijk}^{\mathrm{T}},$

## (42)

$=-{I}_{0}{D}_{ij}\mathrm{diag}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){D}_{ij}^{\mathrm{T}}.$

Now, the Hessian can be written as follows:

## (43)

$\frac{{\partial }^{2}f}{\partial {\sigma }_{t}^{2}}=\sum _{i=1}^{N}\sum _{j=1}^{N}2{I}_{0}^{2}{D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left[{D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\right]}^{\mathrm{T}}\phantom{\rule{0ex}{0ex}}-2{I}_{0}\left({I}_{ij}-{I}_{0}{\mathbf{E}}_{ij}^{\mathrm{T}}{\mathbf{H}}_{ij}\right){D}_{ij}\mathrm{diag}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){D}_{ij}^{\mathrm{T}}.$
Note that ${D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left[{D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\right]}^{\mathrm{T}}$ should not be expanded as ${D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}{D}_{ij}^{\mathrm{T}}$ because it involves a large matrix $\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}$, which is computationally intensive to compute.

## 4.2.4.

#### Discussion

By reusing ${\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}$ and ${D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)$ computed for the Jacobian, the new formulation of Hessian needs:

• $O\left({N}^{2}{M}^{2}\right)$ multiplications for ${D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left[{D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\right]}^{\mathrm{T}}$,

• $O\left(\overline{N}NM\right)$ multiplications for ${D}_{ij}\mathrm{diag}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)$,

• $O\left(\overline{N}{N}^{2}{M}^{2}\right)$ multiplications for ${D}_{ij}\mathrm{diag}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){D}_{ij}^{\mathrm{T}}$,

for each $i$ and $j$. In total, it takes $O\left(\overline{N}{N}^{4}{M}^{2}\right)$ operations to compute ${N}^{2}{M}^{2}$ elements of the Hessian or $O\left(\overline{N}{N}^{2}\right)$ operations per element.

Contrarily, for each element, the previous formulation [Eq. (32)] needs:

• $O\left({\stackrel{‾}{N}}^{2}\right)$ multiplications for ${\mathbf{\gamma }}_{ijx}$,

• $O\left(\stackrel{‾}{N}\right)$ multiplications and the sum for $\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\odot {\mathbf{\beta }}_{ijx}\odot {\mathbf{\beta }}_{ijy}\right)$,

• $O\left({\stackrel{‾}{N}}^{2}\right)$ multiplications for $\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}$,

• $O\left({\stackrel{‾}{N}}^{2}\right)$ multiplications and the sum for $\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\odot {\mathbf{\gamma }}_{ijx}\odot {\mathbf{\gamma }}_{ijy}$,

• $O\left({\stackrel{‾}{N}}^{2}{N}^{2}\right)$ additions for the sums of $i$ and $j$.

In total, it takes $O\left({\stackrel{‾}{N}}^{2}{N}^{2}\right)$ operations to compute a single element of the Hessian.

In summary, the proposed new formulation has the cost of $O\left(\overline{N}{N}^{2}\right)$ operations per element, whereas the previous old formulation has the cost of $O\left({N}^{2}{\stackrel{‾}{N}}^{2}\right)$ operations per element. Table 2 summarizes the discussion above.

## Table 2

Comparison of the new and old formulations for computing the Hessian.

TermsNewOld
${D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left[{D}_{ij}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)\right]}^{\mathrm{T}}$$O\left({N}^{2}{M}^{2}\right)$
${D}_{ij}\mathrm{diag}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)$$O\left(\overline{N}NM\right)$
${D}_{ij}\mathrm{diag}\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){D}_{ij}^{\mathrm{T}}$$O\left(\overline{N}{N}^{2}{M}^{2}\right)$
$\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\odot {\mathbf{\beta }}_{ijx}\odot {\mathbf{\beta }}_{ijy}\right)$$O\left(\stackrel{‾}{N}\right)$
$\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}$$O\left({\overline{N}}^{2}\right)$
$\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right){\left({\mathbf{E}}_{ij}\odot {\mathbf{H}}_{ij}\right)}^{\mathrm{T}}\odot {\mathbf{\gamma }}_{ijx}\odot {\mathbf{\gamma }}_{ijy}$$O\left({\stackrel{‾}{N}}^{2}\right)$
Sums of $i$, $j$$O\left({N}^{2}{\stackrel{‾}{N}}^{2}\right)$
Per element$O\left(\overline{N}{N}^{2}\right)$$O\left({N}^{2}{\overline{N}}^{2}\right)$

## Numerical Simulation

In this section, we report the results obtained by simulations using the proposed method by comparing PD and LB interior point methods, as well as old and new formulations of Jacobian and Hessian.

Since the mathematical model we used to describe the light transport in the forward problem is exactly the same as the model in the previous work,15 we use the same setup as follows. For the 2-D layered medium, the grid size was set to $N=M=24$ with square voxels of size 1 (mm), i.e., the medium is $24\text{\hspace{0.17em}\hspace{0.17em}}\left(\mathrm{mm}\right)×24\text{\hspace{0.17em}\hspace{0.17em}}\left(\mathrm{mm}\right)$. The values of the extinction coefficients are set between 1.05 and 1.55 (${\mathrm{mm}}^{-1}$), and the lower and upper bounds (${\sigma }_{tl}$ and ${\sigma }_{tu}$) are set to be 1.0 and 2.0 (${\mathrm{mm}}^{-1}$), respectively. Values of the initial guess are 1.001 for all elements of ${\mathbf{\sigma }}_{t0}$, as well as ${\mathbf{s}}_{0}$ and ${\mathbf{z}}_{0}$. Parameters used in Algorithm 1 are set as $\sigma =0.5$, $\eta =0.01$, ${\epsilon }_{\mu }=1$, and ${\epsilon }_{\mathrm{TOL}}=0.02$.

## 5.1.

### Estimation Quality

The ground truth and the estimated results are shown in Fig. 3. The matrix plots in the top row of this figure represent five different media (a)–(e) used for the simulation, which were also used in the previous work.15 Note that the medium e is the Shepp–Logan phantom.31 Each voxel is shaded in gray according to the values of the extinction coefficients.

## Fig. 3

Numerical simulation results for five media (a)–(e) in a grid of $24×24$. Darker shades of gray represent larger values of extinction coefficients (more light is absorbed at the voxel). The bars on the side show extinction coefficient values in gray scale. The first row shows ground truth for five different types of media used for the simulation. The following rows show the estimated results for different combinations of LB or PD methods, old or new formulas, and Newton’s or quasi-Newton methods.

## Fig. 4

Numerical simulation results for five media (a)–(e) in a grid of $24×24$. Darker shades of gray represent larger values of extinction coefficients (more light is absorbed at the voxel). The bars on the side show extinction coefficient values in gray scale. The first row shows ground truth for five different types of media used for the simulation. The following rows show the estimated results for different combinations of LB or PD methods, old or new formulas, and Newton’s or quasi-Newton methods. Results of the previous work15 and DOT are shown as baselines in the last three rows.

The following rows show the estimated results of different combinations of LB or PD methods, old or new formulas, and Newton’s or quasi-Newton methods. The proposed method is PD-new-Newton/quasi-Newton; that is, the PD method with Newton’s or quasi-Newton method by using the proposed new formulation. The row LB-old-quasi-Newton corresponds to the previous work15 that uses the LB method with the quasi-Newton method by using the old formulation, and the row LB-old-Newton corresponds to another prior work.12

As we can see, the results of different combinations almost look the same for each of the five media. This observation is also validated by the root-mean-square error (RMSE) shown in Table 3. The RMSE values of all combinations are more or less the same, while some variations appear due to the different update rules between Newton’s and quasi-Newton methods, and different stopping conditions between LB and PD methods.

## Table 3

RMSEs and computation time for the numerical simulations for five different types of media (a)–(e) with grid size of 24×24, for different combinations of LB or PD methods, old or new formulas, and Newton’s or quasi-Newton methods. Each computation time shows the mean and standard deviation of 10 trials, except the combinations of “old-Newton.” Note that RMSE values are exactly the same for 10 trials. Results of DOT methods are shown for comparison.

abcde
RMSELB-new-Newton0.0084220.0126430.0145940.0212460.052511
LB-new-quasi-Newton0.0086460.0124780.0144440.0203750.049811
LB-old-Newton0.0084220.0126430.0145940.0212460.052420
LB-old-quasi-Newton0.0086460.0124780.0144440.0203750.049818
PD-new-Newton0.0097760.0131900.0144900.0212510.055912
PD-new-quasi-Newton0.0097540.0131840.0145020.0212010.056085
PD-old-Newton0.0097760.0131900.0144900.0212510.055912
PD-old-quasi-Newton0.0097540.0131840.0145020.0212010.056084
DOT (GN, Laplace prior)0.0593390.0629840.0781000.0650010.087094
DOT (GN, NOSER prior)0.0520530.0575150.0754780.0591560.086397
DOT (GN, Tikhonov prior)0.0547290.0561960.0731460.0592840.087659
DOT (PD, TV prior)0.0550470.0592190.0818110.0702630.086107
Computation time (s)LB-new-Newton$60.00±4.60$$57.63±1.41$$61.90±2.88$$62.32±1.22$$93.10±2.46$
LB-new-quasi-Newton$18.64±0.90$$17.22±1.03$$25.32±1.21$$20.73±1.10$$32.57±0.58$
LB-old-Newton12610012848133831403721577
LB-old-quasi-Newton$44.86±1.33$$42.58±1.19$$63.75±1.76$$54.05±2.04$$91.17±1.48$
PD-new-Newton$18.73±2.18$$16.52±0.67$$17.40±0.93$$18.28±1.30$$23.14±1.32$
PD-new-quasi-Newton$14.44±0.78$$13.07±0.61$$40.25±1.01$$30.77±0.86$$48.26±1.59$
PD-old-Newton56735418582456637547
PD-old-quasi-Newton$75.67±1.34$$67.18±1.38$$203.42±3.42$$155.69±2.67$$248.05±4.86$
DOT (GN, Laplace prior)$0.34±0.04$$0.41±0.04$$0.40±0.03$$0.40±0.04$$0.40±0.03$
DOT (GN, NOSER prior)$0.52±0.04$$0.52±0.05$$0.53±0.05$$0.50±0.01$$0.52±0.03$
DOT (GN, Tikhonov prior)$0.29±0.01$$0.29±0.02$$0.30±0.03$$0.29±0.01$$0.29±0.02$
DOT (PD, TV prior)$2.67±0.07$$2.68±0.06$$2.67±0.07$$2.64±0.05$$2.66±0.07$

## 5.2.

### Estimation Time

The main goal of this paper is to develop an efficient way to solve the inverse problem. Table 3 shows the computation cost of different combinations. All experiments were performed on a Linux workstation (two Intel Xeon E5-2630 2.4 GHz CPUs, total 16 physical cores, with 256 GB memory). We implemented the method in MATLAB R2017a and did not explicitly use the parallel computation toolbox of MATLAB except the Hessian computation of LB/PD-old-Newton due to its slow computation. Parallel matrix multiplications are, however, automatically performed on MATLAB. Table 3 shows the computation time for different computations in seconds. We report the average and standard deviation of 10 trials, except the cases of LB/PD-old-Newton which show the processing time of a single trial.

For any combination, our proposed new formulation is much faster than the old formulation. The uses of Newton’s method greatly benefit from the efficient Hessian computation and the computation time reduces more than a factor of 100. However, the new formulation does not help to reduce the computation cost of quasi-Newton methods and the reduction is a factor of just 2 or 3. This is due to the fact that the quasi-Newton method needs gradient vectors, and its computation is of the order $NM$ (the number of voxels) and it is not so large in terms of the total computation cost. In contrast, Hessian computation in the Newton’s method is of the size $NM×NM$, which is quite large compared to the gradient. Our new formulation is, therefore, better when the Newton method is used.

With the quasi-Newton method, the PD approach seems to be comparable to the LB method. By comparing rows LB/PD-new-quasi-Newton, LB is faster than PD for denser media (c), (d), and (e). This might be caused by the different ways of approximations by the quasi-Newton method. For the LB method, the gradient is modified by the approximated Hessian. For the PD method, however, the approximated Hessian is used in the matrix to be solved, resulting in an update rule of ${\mathbf{p}}_{\sigma }$ regularized by diagonal elements of $\mathbf{w}$ in Eq. (17). Except the simplest medium (a), the fast combination is PD-new-Newton, which is proposed in this paper. This is due to the fast convergence of Newton’s method compared to the quasi-Newton method, and also the fact that the PD method needs fewer iterations than the LB method. The qualities of results are almost the same as discussed above, then PD-new-Newton is the best when the working memory is enough for storing the Hessian. Otherwise, LB/PD-new-quasi-Newton is better to be used.

## 5.3.

### Comparison

We compare our methods to a standard DOT implemented in the Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software (EIDORS).32,33 in the same setting as the previous work:15 $N=M=24$ medium of size $24\text{\hspace{0.17em}\hspace{0.17em}}\left(\mathrm{mm}\right)×24\text{\hspace{0.17em}\hspace{0.17em}}\left(\mathrm{mm}\right)$ with the five media (a)–(e). For solving DOT by EIDORS, we used $24×24×24=1152$ triangle elements. For boundary conditions, we placed 48 light sources and detectors at the same intervals around the medium. We used some different solvers and priors; Gauss–Newton method34 with Laplace, NOSER,35 and Tikhonov priors, and PD method with total variation prior.

Due to the diffusion approximation of DOT, the results in Fig. 4 for DOT with the Gauss–Newton method are blurred, and those for DOT with PD have a tendency of overestimating the high-coefficient value areas. In contrast, the results of PD-new-Newton (and other combinations in Fig. 3) are clearer and sharper for all combinations. This observation is also validated by the RMSE shown in Table 3. The RMSE values of PD-new-Newton are smaller than the values of DOT for all five media.

The obvious drawback of PD-new-Newton is its high computation cost. It is slower by a factor of 10 compared to DOT with the PD method and a factor of 100 to DOT with the Gauss–Newton method. A large part of the computation cost comes from the computation of Hessian and Jacobian, which depends on the number of paths ${N}_{ij}$. Further reduction of computation cost is left for our future work.

## Conclusion

In this paper, we have proposed a PD approach to optical tomography with a discretized path integral and also efficient formulation for computing Hessian and Jacobian. Numerical simulation examples with 2-D layered media are shown to demonstrate that the proposed method, called PD-new-Newton in the experiments, performed faster than the previous work (LB-old-quasi-Newton) while the estimated extinction coefficients of both methods were comparable. Compared to DOT, the proposed method worked slower but produced better estimation results in terms of RMSE.

## Disclosures

No other conflicts of interest, financial or otherwise, are declared by the authors.

## Acknowledgments

This work was supported in part by the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. JP26280061.

## References

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

2. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. modelling and reconstruction,” Phys. Med. Biol. 42(5), 841–853 (1997).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/42/5/008 Google 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–840 (1997).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/42/5/007 Google Scholar

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

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

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

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

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

9. Y. Mukaigawa, R. Raskar and Y. Yagi, “Analysis of scattering light transport in translucent media,” IPSJ Trans. Comput. Vision Appl. 3, 122–133 (2011). http://dx.doi.org/10.2197/ipsjtcva.3.122 Google Scholar

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

11. Y. Ishii et al., “Scattering tomography by Monte Carlo voting,” in IAPR Int. Conf. on Machine Vision Applications (2013). Google Scholar

12. T. Tamaki et al., “Multiple-scattering optical tomography with layered material,” in 2013 Int. Conf. on Signal-Image Technology & Internet-Based Systems, pp. 93–99, IEEE (2013). http://dx.doi.org/10.1109/SITIS.2013.26 Google Scholar

13. B. Yuan et al., “Layered optical tomography of multiple scattering media with combined constraint optimization,” in 2015 21st Korea-Japan Joint Workshop on Frontiers of Computer Vision (FCV), pp. 1–6, IEEE (2015). http://dx.doi.org/10.1109/FCV.2015.7103735 Google Scholar

14. R. Akashi et al., “Scattering tomography using ellipsoidal mirror,” in 21st Korea-Japan Joint Workshop on Frontiers of Computer Vision (FCV 2015), Mokpo, South Korea, pp. 1–5 (2015). Google Scholar

15. B. Yuan et al., “Optical tomography with discretized path integral,” J. Med. Imaging 2, 033501 (2015).JMEIET0920-5497 http://dx.doi.org/10.1117/1.JMI.2.3.033501 Google Scholar

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

17. B. Yuan et al., “Optical tomography with discretized path integrals: a comparison with log-barrier and primal-dual methods,” in The Korea-Japan joint workshop on Frontiers of Computer Vision (FCV 2016), pp. 378–382 (2016). Google Scholar

18. T. Tamaki et al., “Efficient formulations of optical tomography with discretized path integral,” in The 23th International Workshop on Frontiers of Computer Vision (FCV 2017), pp. 1–5 (2017). Google Scholar

19. G. Zhang et al., “Acceleration of dynamic fluorescence molecular tomography with principal component analysis,” Biomed. Opt. Express 6, 2036 (2015).BOEICL2156-7085 http://dx.doi.org/10.1364/BOE.6.002036 Google Scholar

20. H. M. Huang and I. T. Hsiao, “Accelerating an ordered-subset low-dose X-ray cone beam computed tomography image reconstruction with a power factor and total variation minimization,” PLoS One 11(4), e0153421 (2016).POLNCL1932-6203 http://dx.doi.org/10.1371/journal.pone.0153421 Google Scholar

21. G. Pratx and L. Xing, “GPU computing in medical physics: a review,” Med. Phys. 38(5), 2685–2697 (2011).MPHYA60094-2405 http://dx.doi.org/10.1118/1.3578605 Google Scholar

22. J. Lobera et al., “High performance computing for a 3-D optical diffraction tomographic application in fluid velocimetry,” Opt. Express 23(4), 4021 (2015).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.23.004021 Google Scholar

23. S. Ha et al., “GPU-accelerated forward and back-projections with spatially varying kernels for 3D DIRECT TOF PET reconstruction,” IEEE Trans. Nuclear Sci. 60(1), 166–173 (2013).IETNAE0018-9499 http://dx.doi.org/10.1109/TNS.2012.2233754 Google Scholar

24. F. Xu and K. Mueller, “Accelerating popular tomographic reconstruction algorithms on commodity PC graphics hardware,” IEEE Trans. Nuclear Sci. 52(3), 654–663 (2005).IETNAE0018-9499 http://dx.doi.org/10.1109/TNS.2005.851398 Google Scholar

25. V. G. Nguyen and S. J. Lee, “GPU-accelerated iterative reconstruction from Compton scattered data using a matched pair of conic projector and backprojector,” Comput. Methods Programs Biomed. 131, 27–36 (2016).CMPBEK0169-2607 http://dx.doi.org/10.1016/j.cmpb.2016.04.012 Google Scholar

26. R. Liu, Y. Luo and H. Yu, “GPU-based acceleration for interior tomography,” IEEE Access 2, 841–855 (2014). http://dx.doi.org/10.1109/ACCESS.2014.2349000 Google Scholar

27. M. Schweiger, “GPU-accelerated finite element method for modelling light transport in diffuse optical tomography,” Int. J. Biomed. Imaging 2011, 1–11 (2011). http://dx.doi.org/10.1155/2011/403892 Google Scholar

28. S. Q. Zheng et al., “A distributed multi-GPU system for high speed electron microscopic tomographic reconstruction,” Ultramicroscopy 111(8), 1137–1143 (2011).ULTRD60304-3991 http://dx.doi.org/10.1016/j.ultramic.2011.03.015 Google Scholar

29. A. Borsic, E. A. Attardo and R. J. Halter, “Multi-GPU Jacobian accelerated computing for soft-field tomography,” Physiol. Meas. 33, 1703–1715 (2012).PMEAE30967-3334 http://dx.doi.org/10.1088/0967-3334/33/10/1703 Google Scholar

30. J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed., Springer, New York (2006). Google Scholar

31. L. A. Shepp and B. F. Logan, “The Fourier reconstruction of a head section,” IEEE Trans. Nucl. Sci. 21, 21–43 (1974).IETNAE0018-9499 http://dx.doi.org/10.1109/TNS.1974.6499235 Google Scholar

32. A. Adler and W. R. Lionheart, “EIDORS: towards a community-based extensible software base for EIT,” in 6th Conf. on Biomedical Applications of Electrical Impedance Tomography, London, UK, pp. 1–4 (2005). Google Scholar

33. A. Adler and W. R. Lionheart, “Uses and abuses of EIDORS: an extensible software base for EIT,” Physiol. Meas. 27(5), S25 (2006).PMEAE30967-3334 http://dx.doi.org/10.1088/0967-3334/27/5/S03 Google Scholar

34. A. Adler and R. Guardo, “Electrical impedance tomography: regularized imaging and contrast detection,” IEEE Trans. Med. Imaging 15, 170–179 (1996).ITMID40278-0062 http://dx.doi.org/10.1109/42.491418 Google Scholar

35. M. Cheney et al., “NOSER: an algorithm for solving the inverse conductivity problem,” Int. J. Imaging Syst. Technol. 2(2), 66–75 (1990).IJITEG0899-9457 http://dx.doi.org/10.1002/(ISSN)1098-1098 Google Scholar

## Biography

Bingzhi Yuan received his BE degree in software engineering from 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.

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, Graduate School of Engineering, Hiroshima University, Japan. His current research interests include computer vision, pattern recognition, high-dimensional data analysis, and image processing.

Kazufumi Kaneda received his BE, ME, and DE degrees from Hiroshima University, Japan, in 1982, 1984, and 1991, respectively. He is a professor in the Department of Information Engineering, Graduate School of Engineering, Hiroshima University. 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.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Bingzhi Yuan, Bingzhi Yuan, Toru Tamaki, Toru Tamaki, Bisser Raytchev, Bisser Raytchev, Kazufumi Kaneda, Kazufumi Kaneda, } "Primal-dual approach to optical tomography with discretized path integral with efficient formulations," Journal of Medical Imaging 4(3), 033501 (19 July 2017). https://doi.org/10.1117/1.JMI.4.3.033501 . Submission: Received: 18 April 2017; Accepted: 28 June 2017
Received: 18 April 2017; Accepted: 28 June 2017; Published: 19 July 2017
JOURNAL ARTICLE
10 PAGES

SHARE
KEYWORDS