Deep Learning Can Reverse Photon Migration for
Diffuse Optical Tomography
Jaejun Yooa , Sohail Sabira , Duchang Heob , Kee Hyun Kimb , Abdul Wahaba , Yoonseok Choi c , Seul-I Leed , Eun Young Chaee ,
Hak Hee Kime , Young Min Baeb , Young-wook Choib , Seungryong Choa , and Jong Chul Yea,1
a

Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea.; b Korea Electrotechnology Research Institute (KERI), Ansan 15588, Korea ;
Medical Research Institute, Gangneung Asan Hospital, 38, Bangdong-gil, Sacheon-myeon, Gangneung, Korea; d Asan Biomedical Research Institute, Asan Medical
Center, University of Ulsan College of Medicine, Seoul, Korea; e Department of Radiology and Research Institute of Radiology, Asan Medical Center, University of Ulsan
College of Medicine, Seoul, Korea
c

Can artificial intelligence (AI) learn complicated non-linear physics? Here we propose a novel deep learning approach that learns non-linear
photon scattering physics and obtains accurate 3D distribution of optical anomalies. In contrast to the traditional black-box deep learning
approaches to inverse problems, our deep network learns to invert the Lippmann-Schwinger integral equation which describes the essential
physics of photon migration of diffuse near-infrared (NIR) photons in turbid media. As an example for clinical relevance, we applied the
method to our prototype diffuse optical tomography (DOT). We show that our deep neural network, trained with only simulation data, can
accurately recover the location of anomalies within biomimetic phantoms and live animals without the use of an exogenous contrast agent.
deep learning | diffuse optical tomography | inverse scattering | photon migration

D

RA
FT

eep learning approaches have demonstrated remarkable performance in many computer vision problems, such as
image classification (1) and segmentation (2). With the development of new techniques such as residual network
(3) and rectified linear unit (ReLU) (4), the performance of deep neural networks have continuously improved and
eventually surpassed the human-level-performance in the image classification task (5). Inspired by these successes,
there have been several approaches to apply deep learning methods for bio-medical image reconstruction. Kang et al.
(6) were the first to apply deep convolutional neural network (CNN) in low-dose X-ray computed tomography (CT)
and demonstrated its effectiveness by winning the second place in 2016 AAPM X-ray CT low-dose grand challenge
(7). Using residual learning with U-Net (2) architecture, Han et al. (8) and Jin et al. (9) independently showed that
the globally distributed streaking artifacts from sparse-view X-ray CT can be effectively removed. Deep learning
approaches have been also successfully applied for linearized photo-acoustics (10) and ultrasound imaging (11).
Unlike these imaging applications where the measurement comes from linear operators, there are other imaging
modalities whose imaging physics should be described by complicated non-linear operators. In particular, the diffuse
optical tomography (DOT), which is of our interest in this paper, is notorious due to severely non-linear and ill-posed
operator originated from the diffusive photon migration (12–15). Although near-infrared (NIR) photons can penetrate
several centimeters inside the tissue to allow non-invasive biomedical imaging, the individual photons scatter many
times and migrate along random paths before escaping from or being absorbed by the medium, which makes imaging
task difficult (Fig. 1(a)).
Mathematically, these imaging physics are usually described by partial differential equations (PDE), and the goal
is to recover constitutive parameters of the PDE from the scattered data measured at the boundary. This is called the
inverse scattering problem (16). Classic modalities using the scattered data such as acoustic, electromagnetic or optical
signal falls into this problem. Many dedicated mathematical and computational algorithms for the reconstruction of
location and parameters of anomalies of different geometrical (cavities, cracks, and inclusions) and physical (acoustic,
optical, and elastic) nature have been proposed over the past few decades (17–22). However, most of the classical
techniques are only suited with entire measurements or strong linearization assumptions, which is usually not feasible
in practice.
In a series of papers on diffuse optical tomography (23, 24), electric impedance tomography (25), elastic imaging
problems (26), and optical diffraction tomography (27), our group proposed a novel algorithmic framework which does
not require linearization or iterative updates of Green’s function, yet provides accurate spatial support of inclusions
and their material parameters. This breakthrough comes from a novel interpretation of Lippmann-Schwinger integral
representation of the scattered fields derived in terms of unknown constitutive parameters having jointly sparse
spatial support. Therefore, the support identification problem can be recast as a joint sparse recovery problem
for the unknown densities, invoking various compressed sensing signal recovery algorithms (23–26, 28, 29). Once
the joint supports are identified, the unknown flux and constitute parameters can be estimated by exploiting the
recursive nature of the Lippmann-Schwinger integral equation, which results in a fast and accurate reconstruction.
However, the performance of the method heavily depends on the accurate modeling of the boundary conditions and
acquisition geometry. Therefore, in practical situation where the accurate modeling of the imaging system is difficult,

D

arXiv:1712.00912v1 [cs.CV] 4 Dec 2017

This manuscript was compiled on December 5, 2017

|

December 5, 2017

|

vol. XXX

|

no. XX

|

1–15

Theory

D

RA
FT

the algorithm often produces unsatisfactory results.
In fact, there are a few preliminary works that attempted to solve the inverse scattering problem using machine
learning approaches (30, 31). For example, to obtain a non-linear inverse scattering solution for optical diffraction
tomography, Kamilov et al. (30) proposed so-called beam-propagation method that computes the unknown photon
flux using back-propagation algorithm. The method by Broek and Koch (31) can be also considered an earlier version
of beam-propagation method using neural network to calculate the dynamical scattering of fast electrons. However,
these methods did not consider direct inversion of 3-D distribution of anomalies.
Recently, Ye at al. (32) proposed a novel mathematical framework called deep convolutional framelets to understand
deep learning approaches in inverse problems. The authors found that a deep neural network with encoder-decoder
architecture is closely related to the multi-level decomposition of Hankel matrix constructed from the unknown
signals. In particular, they further showed that Hankel matrix decomposition is equivalent to the convolution framelet
expansion by Yin et al. (33), and training the neural network is actually the process of learning a set of optimal
local-bases for a given non-local bases. Accordingly, imaging application-specific knowledge leads to a better choice of
non-local bases, on which the local bases are learned to maximize the reconstruction performance.
This perspective of deep convolutional framelets gives us an important clue for designing a deep network that
learns non-linear imaging physics. Since the nonlinear mapping of Lippmann-Schwinger type integrator equations
is fundamental in inverse scattering problems as shown in (23–27), this imaging physics should be included in the
design of the network rather than a traditional approach of using the deep network as a black-box. In particular,
our network should be designed to invert the Lippmann-Schwinger equation, but due to the ill-posed nature of
the Lippmann-Schwinger equation, additional requirement for our network is such that the output of the inverse
mapping lies in a low-dimensional manifold. Interestingly, by adding a fully connected layer at the first stage of the
network followed by a CNN with an encoder-decoder structure, this physical intuition is directly mapped to each
layer of convolutional neural network and our neural network can robustly invert the Lippmann-Schwinger equation
to accurately recover the optical anomalies.
As a clinical relevance, we designed a DOT scanner as a part of simultaneous X-ray digital breast tomosynthesis
(DBT) (34) and DOT imaging system, and applied the proposed network architecture as an inversion engine for optical
imaging part. Although the network was trained only using the numerical data generated via diffusion equation,
extensive results using numerical- and real- biomimic phantom as well as in vivo animal experiments substantiate
that the proposed method consistently outperforms the conventional methods.
The paper is organized as follows. We first provide a rigorous derivation of the Lippmann-Schwinger equation from
the diffusion equation. This is followed by a brief introduction to mathematical theory of the deep convolutional
framelet expansion (32) and a discussion on how the inversion of the Lippmann-Schwinger equation can be mapped
to a deep neural network. Then, experimental results follow to confirm the theory. The generality of the proposed
method for other inverse scattering problems will be briefly discussed in Conclusion.

Lippmann-Schwinger Integral Equation for Photon Migration. Diffuse optical tomography uses NIR light to investigate anomalies within the tissue that have different optical absorption and scattering parameters. As these
perturbations of optical parameters are usually due to the hemoglobin concentration, DOT is sensitive to biological
function and disease that are associated with hemoglobin concentration changes. Accordingly, it has been investigated
as a promising non-invasive imaging tool, for example, for functional brain imaging and breast cancer detection
(12–15).
The inverse problem in DOT is to find the locations and geometric features of anomalies Ω1 , · · · , ΩN ⊂ Ω ⊂ R3
with their optical parameters (such as absorption µ1 , · · · , µN and/or diffusion D1 , · · · , DN ), using a finite number of
pairs of optical measurements on a part of its boundary ∂Ω. In particular, our DOT system is mainly interested in
the absorption parameter changes due to the hemoglobin concentration changes represented by

µδ (x) := µ0 χΩ\∪N
(x) +
n=1 Ωn

N
X

µn χΩn (x).

[1]

n=1

whereas the diffusion parameters of the anomalies are indifferent from the background, i.e., Dn = D0 , for all
n = 1, · · · , N . Here, χB represents the characteristic function of a domain B, and the anomalies Ω1 , · · · , ΩN be
The authors declare no conflict of interest.
1

To whom correspondence should be addressed. E-mail: jong.yekaist.ac.kr

2

|

Yoo et al.

compactly supported inside Ω and have boundaries ∂Ω1 , · · · , ∂ΩN , respectively, such that Ω \ ∪N
n=1 Ωn is a simply
connected domain.
If the medium has a characteristic of highly scattering interactions dominating over absorption, which is the case
for the medium such as a cloud milk or biological tissue, then the light propagation can be modeled by the diffusion
equation (12, 18):

∇ · (D0 (x)∇utm (x)) − µδ (x)utm (x) = 0, x ∈ Ω,
t
[2]
utm (x) + ` ∂um (x) = gm (x),
x ∈ ∂Ω,
∂ν
p
where κ0 := µ0 /D0 denotes the diffuse wave-number, utm , for m = 1, · · · , M ∈ N, is the total photon density
corresponding to the photon flux gm at the boundary ∂Ω from the m−th source distribution, ` is the extrapolation
distance, and ∂/∂ν denotes the normal derivative in the direction of outward unit normal ν at ∂Ω.
Let uim be the background photon density for the homogeneous medium Ω (in the absence of any perturbation) and
s
um := utm − uim be the scattered photon flux. Then, the m−th scattered photon flux measurement at the boundary
for the absorption parameter perturbation δµ := µδ − µ0 is given by
usm (x) =Mm [f ](x),

∀x ∈ ∂Ω,

 
f := δµ .

[4]

RA
FT

where the measurement operator Mm is given by
Z
1
Mm [f ](x) := −
G(x, y)utm (y)f (y)dy,
D0 ∪N
Ω
n=1 n

[3]

To facilitate further discussion, we define the forward mapping M by collecting measurements with respect to all
source positions, i.e.,




g := us1 (x) · · · usM (x) = Mf = M1 [f ](x) · · · MM [f ](x) , x ∈ ∂Ωd ,
[5]

D

where ∂Ωd denotes the detector locations at the boundary; see Supplementary Material for further details.
Note that the g in Eq. (5) is a matrix of dimension of (no. of detectors) × (no. of sources). Eq. (3) is the so-called
Lippmann-Schwinger integral equation (derived in Supplementary material) that relates the scattered field to the
perturbation in the absorption parameter. This construction of the multi-static data matrix g was an essential part
of the joint sparse recovery formulations (23–27). However, the measurement operator relating perturbation in the
optical parameter involves the total optical photon density utm (x) which, in turn, also depends on the unknown
perturbation. This makes the inverse problem highly non-linear. Furthermore, due to the dissipative nature of
diffusive wave and smaller number of measurements compared to the number of unknowns, resolving the image
from DOT is a severely ill-posed problem (35). The basic idea of joint sparse recovery formulations (23–27) is then
to decouple the nonlinear inverse problems in two consecutive steps, taking advantage of the fact that the optical
perturbation does not change position during multiple illumination. In this paper, we further extend this idea to the
point that the optical anomalies can be directly recovered from the measured multi-static data matrix g.
Neural network for inverting Lippmann-Schwinger Equation. In this section, we derive a novel neural network to

directly invert the forward operator Eq. (5). Specifically, let f be the 3D distribution of the perturbation δµ in the
region of interest. If f comes from a smoothly varying perturbation, then we can easily expect that the corresponding
Fourier spectrum fˆ(ω) is mostly concentrated in a small number of coefficients, so there exists a spectral domain
function ĥ(ω) such that
fˆ(ω)ĥ(ω) = 0, ∀ω ,
which can be equivalently represented in the discrete spatial domain by
f ~ h = 0,

[6]

where h is called the annihilating filter (36). Eq. (6) implies the existence of a rank deficient block Hankel matrix,
denoted by Hd (f ) ∈ Rn×d , constructed from f and its rank is dependent on the Fourier domain sparsity level as
rigorously shown in (37).
Let the block Hankel matrix Hd (f ) admit the singular value decomposition Hd (f ) = U ΣV > where U = [u1 · · · ur ]
and V = [v1 · · · vr ] denote the left and the right singular vector bases matrices, respectively, Σ ∈ Rr×r is the diagonal
Yoo et al.

|

December 5, 2017

|

vol. XXX

|

no. XX

|

3

Fig. 1. (a) Schematic illustration of photon behavior within the diffusive medium. When photons are illuminated on the tissue, photons can be reflected, scattered or absorbed
in the tissue. The goal of the DOT is to reconstruct the distribution of optical anomalies within tissues from scattered photon measurements. (b) A generic form of a neural
network for inversion of the Lippmann-Schwinger equation. The measured multi-static data matrix is first mapped to the voxel domain by a fully connected layer representing
T , which then generates the framelet coefficient C by the encoder layer filter Ψ. The framelet coefficient C is further filtered using the additional convolution filter H . Finally,
the decoder reconstructs the 3-D distribution of the optical anomalies from the filtered framelet coefficients Ĉ by the decoder layer filter ν(Ψ).

RA
FT

matrix whose diagonal components contain the singular values, and r denotes its rank. If there exist two matrices
pairs (Φ, Φ̃) and (Ψ, Ψ̃) satisfying the conditions
Φ̃Φ> = In×n ,

ΨΨ̃> = PR(V ) ,

[7]

where R(V ) denote the range space of V and PR(V ) represents a projection onto R(V ), then
Hd (f ) = Φ̃Φ> Hd (f )ΨΨ̃> = Φ̃CΨ> ,

[8]

C = Φ> Hd (f )Ψ.

[9]

with the coefficient matrix C given by

D

One of the most important discoveries in (32) is that an encoder-decoder structure convolution layer is emerged from
Eq. (9) and Eq. (8). Precisely, Eq. (8) is equivalent to the following paired encoder-decoder convolution structure:

C = Φ> f ~ Ψ
[10]

f = Φ̃C ~ τ (Ψ̃),
[11]
where the convolutions in Eq. (10) and Eq. (11) corresponds to the single-input multi-output (SIMO) convolution
and multi-input single-output (MISO), respectively:
 
ψ̃1


1 . 
Ψ := ψ 1 · · · ψ r , τ (Ψ̃) :=  ..  ,
[12]
d
ψ̃r
where ψ i denotes the flipped version of the vector ψi , i.e., its indices are reversed (37).
This encoder-decoder representation of the signal is ideally suitable for inverting Lippmann-Schwinger integral
equation for our DOT imaging problem. First, we can simply choose Φ = Φ̃ = I. Then, by defining an inversion
operator T := M−1 for the forward operator in Eq. (5) and substituting f = T g in Eq. (10), the encoder-decoder
structure neural network is re-written as
 
C = (T g) ~ Ψ,
f = Ĉ ~ τ (Ψ̃)
[13]
where Ĉ could be further processed coefficient from C to remove noise, i.e.
Ĉ = C ~ H,
4

|

Yoo et al.

Fig. 2. Multi-channel DOT system used in our experiments (38, 39). The light source has three fiber pigtailed laser diode modules with 785 nm, 808nm and 850 nm. 70
MHz RF signal is simultaneously applied to these light sources using bias-T, RF splitter and RF AMP. Two optical switches are used to deliver light to 64 specific positions
in the source probe. During optical switching time, one-tone modulation light photons reach 40 detection fiber ends after passing an optical phantom and are detected
simultaneously by 40 avalanche photodiodes (APD) installed in the home-made signal processing card.

D

RA
FT

for some filter H. Finally, by choosing appropriate number of the filter channel r, we can control the dimension of the
reconstruction image manifold.
The corresponding three-layer network structure is illustrated in Fig. 1(b). Here, the network consists of a
single fully connected layer that represents T , two paired 3D-convolutional layers with filters Ψ and τ (Ψ̃), and the
intermediate 3D-convolutional layer for additional filtering. To enable learning from the training data, we used the
hyperbolic tangent function (tanh) as an activation function for the fully connected layer and two convolutional layers
(C1 and C2), whereas the last convolutional layer (C3) was combined with rectified linear unit (ReLU) to ensure the
positive value for the optical property distribution. Therefore, the network is designed such that the inversion of
the Lippmann-Schwinger equation is performed to lead the resulting optical parameters in low-dimensional manifold
structure.
Compared to the conventional model-based approaches, the proposed deep network has many advantages. First,
the inversion of the Lippmann-Swinger equation is fully data-driven such that we do not need any explicit modeling
of the acquisition system and boundary conditions. Second, unlike the explicit regularization in the model-based
approaches, the low-dimensional manifold structure of the desired optical distribution is embedded as convolutional
layers that are also learned from the data. This is why the proposed neural network is felicitous to provide a robust
reconstruction. This will be substantiated through extensive experiments.

Methods
DOT Hardware System. Fig. 2 shows the schematics of the frequency domain diffuse optical tomography (DOT)

system, which we used in this study. The DOT system has been developed at the Korea Electrotechnology Research
Institute (KERI) to improve diagnostic accuracy of the digital breast tomosynthesis (DBT) system by combining the
DOT system with the DBT system for joint breast cancer diagnosis(38, 39). After its development, the system has
been installed and been under clinical test at Asan Medical Center (AMC) (40). The DOT system briefly consists of
four parts: light source, optical detector, optical probe, and data acquisition and controller. The light source has
three fiber pigtailed laser diode modules with 785 nm, 808nm and 850 nm. 70 MHz RF signal is simultaneously
applied to these light sources using bias-T, RF splitter and RF AMP. Two optical switches are used to deliver light to
64 specific positions in the source probe. During optical switching time, one-tone modulation light photons reach 40
detection fiber ends after passing an optical phantom and are detected simultaneously by 40 avalanche photodiodes
(APD) installed in the home-made signal processing card. The DOT system uses an In-phase(I) and quadrature(Q)
demodulator to get amplitude and phase of the signal in the signal processing card. The 40 IQ signal pairs are
simultaneously acquired using data acquisition boards.
Biomimic phantoms. Next, to analyze the performance of the proposed approach in controlled real experiments, two

biomimic phantoms with known inhomogeneity locations are created (see Fig. 4). The first phantom is made of
polypropylene containing a vertically oriented cylindrical cavity that has 20 mm diameter and 15 mm height. The
Yoo et al.

|

December 5, 2017

|

vol. XXX

|

no. XX

|

5

cavity is filled with the acetyl inclusion which has different optical properties to the background. For the second
phantom, we used a custom-made open-top acrylic chamber (175mm x 120mm x 40mm) and three different sized
knots (5mm, 10mm, 20mm diameter) for the mimicry of a tumor-like vascular structure. The knots were made using
thin polymer tube (I.D 0.40mm, O.D 0.8mm diameter) and were filled with the rodent blood that was originated from
the abdominal aorta of Sprague-Dawley rat that was under 1 to 2% isoflurane inhalation anesthesia. The chamber was
filled with the completely melted pig lard and the medium was coagulated at room temperature for the imaging scan.
Animal imaging. The mouse colon cancer cell line MC38 were obtained from Scripps Korea Antibody Institute(Chuncheon,

Korea) and the cell line was cultivated in Dulbecco’s modified Eagle’s medium (DMEM, GIBCO, NY, US) supplemented with 10% Fetal bovine serum(FBS, GIBCO) and 1x Antibiotic-Antimycotic(GIBCO). For the tumor-bearing
mice, 5x106 cells were injected subcutaneously into the right flank region of C57BL/6 mice aging 7-9 weeks (Orient
Bio, Seongnam, Korea). Animal hairs were removed through trimming and waxing. Anesthesia was applied during
the imaging scanning with an intramuscular injection of Zoletil and Rumpun (4:1 ratio) in normal saline solution.
Mice were placed inside of the custom-made 80mm x 80mm x 30mm open-top acrylic chamber that had a semicircle
hole-structure on the one side of the chamber for the relaxed breathing. A gap between the semicircle structure
and the head was sealed with the clay. The chamber was filled with the water/milk mixture as 1000:50 ratios. All
experiments associated with this study were approved by Institutional Animal Care and Use Committees of Asan
Medical Center (IACUC no. 2017-12-198).
Data preprocessing. To determine the maximally usable source-detector distance, we measured signal magnitude

RA
FT

according to the source-detector distances. We observed that the signals were above the noise floor when the separation
distance between the source and the detector ρ was less than 51 mm (ρ < 51 mm). Therefore, instead of using
measurements at all source-detector pairs, we only used the pairs having source and detector less than 51 mm apart.
This step not only enhanced the signal-to-noise ratio (SNR) of the data but also largely reduced the number of
parameters to be learned in the fully connected layer. For example, in the source-detector configuration of the
biomimic phantom, the number of source-detector pairs reduced from 2560 to 466. This decreased the number of the
parameters to train from 137,625,600 to 25,105,920, which is an order difference. To match the scale and bias of the
signal amplitude between the simulation and the real data, we multiplied appropriate calibration factor to the real
measurement to match the signal envelope from the simulation data. We centered the data cloud on the origin with
the maximum width of one by subtracting the mean across every individual data and dividing it with its maximum
value. To deal with the unbalanced distribution of nonzero values in the 3D label image and to prevent the proposed
network learning a trivial mapping (rendering all zero values), we weighted the non-zero values by multiplying a
constant according to the ratio of the total voxel numbers over the non-zero voxels.

D

Neural network training. In order to test the robustness of the deep network in real experiments and to obtain a large

database in an efficient manner, the training data were generated by solving Eq. (2) using finite element method
(FEM) based solver NIRFAST (see, e.g., (41, 42)). The finite element meshes were constructed according to the
specifications of the phantom used in each experiment (see Table 1). We generated 1500 numbers of data by randomly
adding up to three spherical heterogeneities of different sizes (having radius between 2 mm to 13 mm) and optical
properties in the homogeneous background. The optical parameters of the heterogeneities were constrained to lie in a
biologically relevant range, i.e., two to five times bigger than the background values. The source-detector configuration
of the data is set to match that of real experimental data displayed in Table 1. To make the label in a matrix form,
FEM mesh is converted to the matrix of voxels by using triangulation-based nearest neighbor interpolation with
an in-built MATLAB griddata function. The number of voxels per each dimension used for each experiment can
be found in Table 1. To train the network for different sizes of phantoms and source-detector configurations, we
generated different sets of training data and changed the input and output sizes of the network, accordingly. The
specifications of the network architecture are provided in Table 2. Note that the overall structure of the network
remained the same except the specific parameters.
The input of the neural network is the multi-static data matrix of pre-processed measurements. To perform
convolution and to match its dimension with the final output of 3D image, the output of the fully connected layer is set
to the size of the discretized dimension for each phantom. All the convolutional layers were preceded by appropriate
zero padding to preserve the size. We used tanh as a non-linear activation function for every layer except the last
convolutional layer where ReLU is used. In the network structure for polypropylene phantom, for example, the first
convolutional layer convolves 16 filters of 3 × 3 × 3 with stride 1 followed by tanh. The second convolutional layer
again convolves the filter of 3 × 3 × 3 with stride 1 followed by ReLU (Table 2).
We used the mean squared error (MSE) as a loss function and the network was implemented using Keras library
(43). Weights for all the convolutional layers were initialized using Xavier initialization. We divided the generated
6

|

Yoo et al.

Table 1. Specification of FEM mesh for test data generation.

# of sources
Polyprophylene phantom
Biomimic phantom
Mouse(normal)
Mouse (with tumor)

# of detectors

4 × 16
4 × 16
7×4
7×7

nodes
20,609
53,760
12,288
12,288

5×8
5×8
7×4
7×7

Optical parameters
Background (mm−1 )

FEM mesh
elements
86,284
291,870
63,426
63,426

µa

µ0s

0.003
0.002
0.0041
0.0045

0.5
1
0.4503
0.3452

# of voxels per xyz dimensions
(resolution= 2.5 mm)

32 × 64 × 20
48 × 70 × 16
32 × 32 × 12
32 × 32 × 12

Table 2. Network architecture specifications. Here, N M is the number of filtered measurement pairs (Polypropylene: N M =
538, Biomimic: N M = 466, Mouse (normal): N M = 470, Mouse (with tumor): N M = 1533).
Type

Network input
Fully connected
dropout
3D convolution
3D convolution
3D convolution

patch size
/stride
3 × 3 × 3/1
3 × 3 × 3/1
3 × 3 × 3/1

Polypropylene
output
size

depth

1 × NM
32 × 64 × 20 × 1
-

32 × 64 × 20 × 16
32 × 64 × 20 × 16
32 × 64 × 20 × 1

16
16
1

patch size
/stride
3 × 3 × 3/1
3 × 3 × 3/1
3 × 3 × 3/1

Biomimic
output
size

depth

1 × NM
48 × 70 × 16 × 1

64
64
1

-

48 × 70 × 16 × 64
48 × 70 × 16 × 64
48 × 70 × 16 × 1

patch size
/stride
3 × 3 × 3/1
3 × 3 × 3/1
3 × 3 × 3/1

Animal
output
size

depth

1 × NM
32 × 32 × 12 × 2

128
128
1

-

32 × 32 × 12 × 128
32 × 32 × 12 × 128
32 × 32 × 12 × 1

D

RA
FT

data into 1000 training and 500 validation data sets. For training, we used the batch size of 64 and Adam optimizer
(44) with the default parameters as mentioned in the original paper, i.e., we used learning rate= 0.0001, β1 = 0.9, and
β2 = 0.999. Training runs for up to 120 epochs with early stopping if the validation loss has not improved in the
last 10 epochs. To prevent overfitting, we added a zero-centered Gaussian noise with standard deviation σ = 0.2
and applied dropout on the fully connected layer with probability p = 0.7. We used a GTX 1080 graphic processor
and i7-6700 CPU (3.40 GHz). The network took about 380 seconds for training. Since our network only used the
generated simulation data for training, it suffered from the noise which cannot be observed from the synthetic data
and overfitting due to lack of complexity. However, by carefully matching the signal envelops of the simulation data
to those of the real data and tuning the parameters of the modules such as dropout probability and the standard
deviation σ of the Gaussian noise, we could achieve the current network architecture which performs well in various
experimental situations. We have not imposed any other augmentation such as shifting and tilting since our input
data is not in the image domain but in the measurement domain which is unreasonable to apply such techniques.
Every 3D visualization of the results is done by using ParaView (45).
Baseline algorithm for comparison. The baseline iterative algorithm was developed using the Rytov iterative method

similar to (21). Specifically, after estimating the bulk optical properties as an initial guess, the Lippmann-Schwinger
equation is linearized using Rytov approximation, and the associated penalized least squares optimization problem
with l2 penalty was solved. The optical properties are updated until the algorithm converges, and the background
Green’s functions are updated accordingly so that the penalized least squares solutions are obtained using newly
updated Green’s function at each iteration. We set the convergence criterion if the reconstructed optical parameter
at current iteration has not improved in the last two iterations. Unless an initial guess is not closer, the algorithm
generally converged in six to ten iterations and each iteration took approximately 40 seconds, which makes total
working time less than 10 minutes.

Results
To validate our deep learning approach in a quantitative manner, the reconstruction experiments from numerical
phantom were first performed, and some representative results are shown in Fig. 3. Here, we used the network trained
using the biomimic phantom geometry (see Table 2). The ground truth images are visualized with binary values
to show the location of virtual anomalies clearly. To demonstrate the robustness of the algorithm under various
situations, we randomly chose the data with different anomaly sizes and locations, with distinct z-locations. Our
proposed network successfully resolved the anomaly locations.
Next, to analyze the performance of the proposed approach in controlled real experiments, two biomimic phantoms
with known inhomogeneity locations are created (see Fig. 4). The first phantom is made of polypropylene containing a
vertically oriented cylindrical cavity that has 20 mm diameter and 15 mm height. The cavity is filled with the acetyl
inclusion which has different optical properties to the background. The schematic illustrations of the phantoms are
Yoo et al.

|

December 5, 2017

|

vol. XXX

|

no. XX

|

7

Fig. 3. The reconstruction results from numerical phantom. Here, we used the network trained using the biomimic phantom geometry (see Table 2). The ground truth images
are visualized with binary values to show the location of virtual anomalies clearly. Our proposed network successfully resolved the anomaly locations.

D

RA
FT

shown with their specifications in Fig. 4. Then, we obtained the measurement data using our multi-channel system.
The reconstructed 3D images from the conventional iterative method and our proposed network are compared. For
the phantom with a single inclusion, both the conventional iterative method and our network accurately reconstructed
the location of optical anomalies (Fig. 4 (top row)). The contrast is also clearly seen in the DBT image. The second
phantom is made of lards to mimic a condition similar to the human breast. We inserted tubes inside the phantom
and made three knots of different sizes. After solidifying the lards, the tubes were filled with blood to mimic a
tumor-like vascular structure. In the DBT image, the locations of the knotted tubes are barely seen, as pointed out
by the black arrows. In this phantom case, the reconstructed image using the iterative method exhibited globally
distributed artifacts, and suffered from artifacts even after thresholding out the values smaller than the 60% of the
max value (Fig. 4 (bottom row)). Since the artifacts had higher values, which dominate the signal from the inclusions,
we manually stopped the algorithm at iteration three to prevent the smoothing regularization from erasing the real
signal. In contrast, the reconstructed image using the proposed network finds three inclusions without those artifacts
showing a better contrast than DBT image. Furthermore, it is remarkable that our reconstruction has a good contrast
so that the knots are easily seen, whereas it was hard to distinguish the virtual tumors from the background using
the DBT image.
Finally, we performed in-vivo animal experiments using a mouse (Fig. 5 and Fig. 6). In order to get the scattered
photon density measurements, we recorded the data with and without the mouse placed in the chamber, which is
filled with the water/milk mixture of 1000 : 50 ratio. Optical scattering data from animal experiments were collected
using the single-channel system. Fig. 5 and Fig. 6 shows the reconstructed images of the mouse with and without
tumor. Both the conventional and the proposed methods recovered high µa values at the chest area of the mouse.
However, the iterative reconstruction finds a big chunk of high µa values around the left thigh of the mouse, which is
unlikely with the normal mouse (Fig. 5). In contrast, our proposed network shows a high µa along the spine of the
mouse where the artery and organs are located. Furthermore, in the mouse with tumor case, our network finds a high
µa values around the right thigh, where the tumor is located (Fig. 6). The lateral view of our reconstructed images
also match with the actual position of the mouse, whose head and body are held a little above the bottom plane due
to the experiment setup.
Compared to the results of the conventional iterative method, our network showed a robust performance over the
various examples. While the iterative reconstruction algorithm often imposed high intensities on spurious locations,
such as the bottom plane of the phantom (Fig. 4 bottom row) and the big clutter at the left thigh of the normal
mouse (Fig. 5 top row), our network found accurate positions with high values only at the locations where inclusions
are likely to exist. Unlike the conventional method, which requires a parameter tuning for every individual case, it
can infer from the measured data without additional pre- and post-processing techniques, which is desirable in the
practical applications.
Note that our network had not seen any real data during the training nor the validation process. However, it
successfully finds the inversion by learning only from the simulation data. Furthermore, even though we trained the
network with sparse examples only, our network successfully finds both sparse (phantom, Fig. 4) and extended targets
(mouse, Fig. 5 and Fig. 6) without any help of regularizers.
To further show that the proposed architecture is near-optimal, we performed ablation studies by changing or
removing the components of the proposed architecture. Since our output f is a 3D distribution, the network needs
to find a set of 3D filters Ψ̄ and τ (Ψ̃). We observed that the network with 3D-convolution showed better z−axis
8

|

Yoo et al.

RA
FT

D

Fig. 4. The reconstructed images of the phantoms. A phantom with a single acetyl cylindrical inclusion (top row) and a biomimic phantom with three blood inclusions (middle
and bottom rows) were used. The first phantom is made of polypropylene containing a vertically oriented cylindrical cavity that has 20 mm diameter and 15 mm height. The
cavity is filled with the acetyl inclusion which has different optical properties to the background. For this experiment, both the conventional iterative method and our network
accurately reconstructed the location of optical anomalies. The contrast is also clearly seen in the DBT image. The second phantom is made of lards to mimic a condition
similar to the human breast. We inserted tubes inside the phantom and made three knots of different sizes. After solidifying the lards, the tubes were filled with blood to mimic
a tumor-like vascular structure. In the DBT image, the locations of the knotted tubes are barely seen, as pointed out by the black arrows. In this case, the reconstructed
image using the iterative method suffered from artifacts at the bottom plane near the detector array. In contrast, the reconstructed image using the proposed network finds
three inclusions without those artifacts showing a better contrast than DBT image.

resolution compared to the one using 2D convolution (Fig. 7), which is consistent with our theoretical prediction.
One may suspect that the performance of the network has originated solely from the first layer since over 98% of
the trainable parameters are from the fully connected layer. To address this concern, we tested the network with
and without convolutional layers after the fully connected layer. We observed that the performance of our network
deteriorated severely and it failed to train without the consecutive convolution layers. At least a single convolution
layer with a single filter were needed to recover the accurate location of the optical anomalies (Fig. 7). However,
paired encoder-decoder filter in the proposed network is better than just using a single convolution layer, which
coincides with our theoretical prediction.
Finally, we observed that the network training become inaccurate if the activation functions do not match with the
physical intuition (Fig. 8). Since the first fully connected layer should learn the non-linear inversion of LippmannSchwinger equation T and the consecutive 3D-convolution layers are supposed to learn the local bases Ψ̄ and their
adjoints τ (Ψ̃), it appeared unreasonable to put non-negativity by imposing a ReLU activation function at each layer.
Therefore, we used tanh as an activation function for every layer except the last convolutional layer outputting the
3D distribution of absorption coefficients which are positive values. This gave the most stable and accurate training.
On the other hand, if we change the activation functions of the fully connected layer and intermediate layers from
tanh to ReLU, we found the network is hardly trained or the performance is degraded, which runs counter to the
common sense in computer vision area that ReLU gives a robust training. Interestingly, if we only have ReLU for
the 3D-convolutional layers, we could observe that the training is very unstable (Fig. 8). Only after the activation
functions of both fully connected layer and the following 3D-convolution layer have tanh activation, the network
Yoo et al.

|

December 5, 2017

|

vol. XXX

|

no. XX

|

9

D

RA
FT

Fig. 5. The reconstructed images of in vivo animal experiments. Mouse without any tumor is shown. In order to get the scattered photon density measurements, we recorded
the data with and without the mouse placed in the tank, which is filled with the water/milk mixture of 1000 : 50 ratio. Both the conventional and the proposed methods
recovered high µa values at the chest area of the mouse. However, the iterative reconstruction finds a big chunk of high µa values around the left thigh of the mouse, which
is unlikely with the normal mouse. In contrast, our proposed network shows a high µa along the spine of the mouse where the artery and organs are located.

Fig. 6. The reconstructed images of in vivo animal experiments. Mouse with tumor on the right thigh is shown. Comparison between the mouse with and without tumor 3D
visualizations are displayed. In the 3D images, the values are thresholded above 60% of the max value for a clear visualization. Our network finds a high µa values around
the right thigh, where the tumor is located (Fig. 6).

starts to learn and shows a stable training (Fig. 8).

Conclusion
In this paper, we proposed a deep learning approach to solve the inverse scattering problem of diffuse optical
tomography (DOT). Unlike the conventional deep learning approach,which tries to denoise or remove the artifacts
from image to image using a black-box approach for neural network, our network was designed based on LippmannSchwinger equation to learn the complicated non-linear physics of the inverse scattering problem. Even though
our network was only trained on the generated numerical data, we showed that the learned mapping is general
over the real experimental data. The simulation and real experimental results showed that the proposed method
outperforms the conventional method and accurately reconstructs the anomalies without iterative procedure or linear
approximation. By using our deep learning framework, the non-linear inverse problem of DOT can be solved in
end-to-end fashion and new data can be efficiently processed in a few hundreds of milliseconds, so it would be useful

10

|

Yoo et al.

RA
FT

Fig. 7. Network ablation study results. The reconstructed images from the networks with 2D- and 3D-convolution are compared (the 2nd and 3rd columns). To show the
necessity of the convolutional layer, the image reconstructed using the network with a single 3D-convolution layer of a single filter is shown (last column). Meanwhile, the
network with only a fully connected layer failed to train (results not shown).

D

Fig. 8. The role of activation functions. The learning curves of the proposed network using tanh as an activation function for every layer except the last ReLU are shown
(the 1st column). The learning curves of the networks with changed activation functions are shown. If we change the activation functions of the fully connected layer and
intermediate layers from tanh to ReLU, we found the network is hardly trained or the performance is degraded (the 2nd column). If we only have ReLU for the 3D-convolutional
layers, we could observe that the training is very unstable (the third column).

for dynamic imaging applications. Moreover, our new design principle based on deep convolutional framelets is a
quite general framework to incorporate physics in network design, so we believe that the proposed network can be
used for other inverse scattering problems from complicated physics.
1. Krizhevsky A, Sutskever I, Hinton GE (2012) Imagenet classification with deep convolutional neural networks in Advances in neural information processing systems. pp. 1097–1105.
2. Ronneberger O, Fischer P, Brox T (2015) U-net: Convolutional networks for biomedical image segmentation in International Conference on Medical Image Computing and Computer-Assisted
Intervention. (Springer), pp. 234–241.
3. He K, Zhang X, Ren S, Sun J (2016) Deep residual learning for image recognition in Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 770–778.
4. Nair V, Hinton GE (2010) Rectified linear units improve restricted boltzmann machines in Proceedings of the 27th international conference on machine learning (ICML-10). pp. 807–814.
5. Russakovsky O, et al. (2015) ImageNet Large Scale Visual Recognition Challenge. International Journal of Computer Vision (IJCV) 115(3):211–252.
6. Kang E, Min J, Ye JC (2017) A deep convolutional neural network using directional wavelets for low-dose x-ray ct reconstruction. Medical Physics 44(10).
7. McCollough CH, et al. (2017) Low-dose ct for the detection and classification of metastatic liver lesions: Results of the 2016 low dose ct grand challenge. Medical Physics 44(10).
8. Han Y, Yoo J, Ye JC (2016) Deep residual learning for compressed sensing ct reconstruction via persistent homology analysis. arXiv preprint arXiv:1611.06391.
9. Jin KH, McCann MT, Froustey E, Unser M (2017) Deep convolutional neural network for inverse problems in imaging. IEEE Transactions on Image Processing 26(9):4509–4522.
10. Antholzer S, Haltmeier M, Schwab J (2017) Deep Learning for Photoacoustic Tomography from Sparse Data. arXiv preprint arXiv:1704.04587.
11. Yoon YH, Ye JC (2017) Deep learning for accelerated ultrasound imaging. arXiv preprint arXiv:1710.10006.
12. Yodh A, Chance B (1995) Spectroscopy and imaging with diffusing light. Physics Today 48:34–40.
13. Ntziachristos V, Yodh A, Schnall M, Chance B (2000) Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement. Proceedings of the National Academy of
Sciences 97(6):2767–2772.
14. Strangman G, Boas DA, Sutton JP (2002) Non-invasive neuroimaging using near-infrared light. Biological psychiatry 52(7):679–693.
15. Boas DA, et al. (2001) Imaging the body with diffuse optical tomography. IEEE signal processing magazine 18(6):57–75.
16. Colton D, Kress R (2013) Integral equation methods in scattering theory. (SIAM).
17. Ntziachristos V, Weissleder R (2001) Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation. Optics letters 26(12):893–895.
18. Markel VA, Schotland JC (2001) Inverse problem in optical diffusion tomography. I. Fourier–Laplace inversion formulas. JOSA A 18(6):1336–1347.
19. Ammari H, et al. (2015) Mathematical methods in elasticity imaging. (Princeton University Press).
20. Jiang H, Pogue BW, Patterson MS, Paulsen KD, Osterberg UL (1996) Optical image reconstruction using frequency-domain data: simulations and experiments. JOSA A 13(2):253–266.
21. Yao Y, Wang Y, Pei Y, Zhu W, Barbour RL (1997) Frequency-domain optical imaging of absorption and scattering distributions by a Born iterative method. JOSA A 14(1):325–342.
22. Ye JC, Webb KJ, Bouman CA, Millane RP (1999) Optical diffusion tomography by iterative-coordinate-descent optimization in a Bayesian framework. JOSA A 16(10):2400–2412.

Yoo et al.

|

December 5, 2017

|

vol. XXX

|

no. XX

|

11

23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.

D

RA
FT

42.
43.
44.
45.
46.

Lee O, Kim JM, Bresler Y, Ye JC (2011) Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity. IEEE transactions on medical imaging 30(5):1129–1142.
Lee O, Ye JC (2013) Joint sparsity-driven non-iterative simultaneous reconstruction of absorption and scattering in diffuse optical tomography. Optics express 21(22):26589–26604.
Lee OK, Kang H, Ye JC, Lim M (2015) A non-iterative method for the electrical impedance tomography based on joint sparse recovery. Inverse Problems 31(7):075002.
Yoo J, Jung Y, Lim M, Ye JC, Wahab A (2017) A joint sparse recovery framework for accurate reconstruction of inclusions in elastic media. SIAM Journal on Imaging Sciences 10(3):1104–1138.
Lim J, et al. (2017) Beyond Born-Rytov limit for super-resolution optical diffraction tomography. Optics express (in press).
Ye JC, Lee SY, Bresler Y (2008) Exact reconstruction formula for diffuse optical tomography using simultaneous sparse representation in Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008.
5th IEEE International Symposium on. (IEEE), pp. 1621–1624.
Ye JC, Lee SY (2008) Non-iterative exact inverse scattering using simultaneous orthogonal matching pursuit (S-OMP) in Acoustics, Speech and Signal Processing, 2008. ICASSP 2008. IEEE
International Conference on. (IEEE), pp. 2457–2460.
Kamilov US, et al. (2015) Learning approach to optical tomography. Optica 2(6):517–522.
Van den Broek W, Koch CT (2012) Method for retrieval of the three-dimensional object potential by inversion of dynamical electron scattering. Physical review letters 109(24):245502.
Ye JC, Han YS (2017) Deep convolutional framelets: A general deep learning for inverse problems. arXiv preprint arXiv:1707.00372.
Yin R, Gao T, Lu YM, Daubechies I (2017) A tale of two bases: Local-nonlocal regularization on image patches with convolution framelets. SIAM Journal on Imaging Sciences 10(2):711–750.
Poplack SP, Tosteson TD, Kogel CA, Nagy HM (2007) Digital breast tomosynthesis: initial experience in 98 women with abnormal digital screening mammography. American Journal of Roentgenology
189(3):616–623.
Arridge SR (1999) Optical tomography in medical imaging. Inverse problems 15(2):R41.
Vetterli M, Marziliano P, Blu T (2002) Sampling signals with finite rate of innovation. IEEE Transactions on Signal Processing 50(6):1417–1428.
Ye JC, Kim JM, Jin KH, Lee K (2017) Compressive sampling using annihilating filter-based low-rank interpolation. IEEE Transactions on Information Theory 63(2):777–801.
Heo D, et al. (2017) Dual modality breast imaging system: Combination digital breast tomosynthesis and diffuse optical tomography. SPIE Photonics west pp. 10059–67.
Heo D, et al. (2017) Dual modality breast imaging system: Combination digital breast tomosynthesis and diffuse optical tomography. Korean Physical Society Spring Meeting.
Choi YW, Park Hs, Kim Ys, Kim HL, Choi JG (2012) Effect of Acquisition Parameters on Digital Breast Tomosynthesis: Total Angular Range and Number of Projection Views. Korean Physical
Society 61(11):1877–1883.
Dehghani H, et al. (2009) Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction. International Journal for Numerical Methods in Biomedical
Engineering 25(6):711–732.
Paulsen KD, Jiang H (1995) Spatially varying optical property reconstruction using a finite element diffusion equation approximation. Medical Physics 22(6):691–701.
Chollet F (2015) Keras (https://github.com/fchollet/keras).
Kingma D, Ba J (2014) Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
Ahrens J, Geveci B, Law C (2005) Paraview: An end-user tool for large data visualization. The Visualization Handbook 717.
Pogue BW, Paulsen KD, Abele C, Kaufman H (2000) Calibration of near-infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantitation in neonatal head-simulating
phantoms. Journal of biomedical optics 5(2):185–193.

12

|

Yoo et al.

Supplementary Material
A. Derivation of the Lippmann-Schwinger Equation for DOT
In this appendix, we derive the general form of the Lippmann-Schwinger equation if both absorption and scattering
changes are present. Although the derivation of the Lippmann-Schwinger equation is already available in the
literature(20–22), the derivation is mainly based on physical intuition rather than mathematical rigor. Therefore, we
provide here a more rigorous mathematical derivation. Let the optical parameter distribution within the domain of
interest Ω be described by
µδ (x) := µ0 χΩ\∪N
(x) +
n=1 Ωn

N
X

µn χΩn (x) and Dδ (x) := D0 χΩ\∪N
(x) +
n=1 Ωn

n=1

N
X

Dn χΩn (x),

[14]

n=1

where χB represents the characteristic function of a domain B. We define the perturbations in optical parameters by
δµ := µδ − µ0

and δD := Dδ − D0 .

x ∈ Ω \ ∪N
n=1 Ωn ,
x ∈ Ωn , n = 1, · · · , N,
[16]

(x), x ∈ ∂Ωn , ∀n = 1, · · · , N,

RA
FT

Then, the total photon density satisfies the transmission problem

t
2 t

∆um (x) − κ0 um (x) = 0,



∇ · (Dn (x)∇utm (x)) − µn (x)utm (x) = 0,

∂utm
∂utm
t
t
u
(x)
=
u
(x)
and
D
(x)
=
D
0
n

m
m
−
+

∂ν −
∂ν


t


ut (x) + ` ∂um (x) = gm (x),
m
∂ν

[15]

+

x ∈ ∂Ω,

D

where w|± (x) := lim→0+ w(x ± ν(x)) for any function w and x ∈ ∂Ωn .
As for mathematical formalism to derive the Lippmann-Schwinger equation for the scattered photon density
usm := utm − uim , the Robin Green function (G), and the background photon density (uim ) for the homogeneous
medium Ω (in the absence of any anomaly) are required. The functions G and uim are, respectively, the solutions to
the boundary value problems

∆x G(x, y) − κ20 G(x, y) = −δy (x), x ∈ Ω,
[17]
∂
G(x, y) + `
G(x, y) = 0,
x ∈ ∂Ω,
∂ν(x)
and


∆uim (x) − κ20 uim (x) = 0,
x ∈ Ω,
i
∂u
uim (x) + ` m (x) = gm (x), x ∈ ∂Ω.
∂ν

[18]

Notice that, by an appropriate use of the Green’s theorem and the boundary condition in Eq. (16),


Z
Z
1
1
∂ut
uim (x) =
G(x, y)gm (y)dσ(y) =
G(x, y) utm (y) + ` m (y) dσ(y), ∀ x ∈ Ω,
` ∂Ω
` ∂Ω
∂ν

[19]

where dσ is the infinitesimal differential element on ∂Ω. Moreover, thanks to the boundary condition in Eq. (17),

Z 
1
∂G
G(x, y) + `
(x, y) utm (y)dσ(y) = 0, ∀ x ∈ Ω.
[20]
` ∂Ω
∂ν(y)
Therefore, by combining Eq. (19) and Eq. (20), it can be easily seen that

Z 
∂utm
∂
t
i
um (x) =
G(x, y)
(y) −
G(x, y)um (y) dσ(y),
∂ν
∂ν(y)
∂Ω

∀ x ∈ Ω.

From an application of the Green’s theorem on Ω \ ∪N
n=1 Ωn , one arrives at
uim (x)

Z
=

Yoo et al.

Ω\∪N
n=1 Ωn

G(x, y)∆utm (y)

−

∆y G(x, y)utm (y)



dy +

N Z
X
n=1

∂Ωn



∂ut
G(x, y) m
∂ν
|

∂
(y) −
G(x, y)utm
∂ν(y)
+

December 5, 2017

|

vol. XXX

|

no. XX


(y) dσ(y).
+

|

13

Thanks to Eq. (16), one can see that
uim (x) = −

Z
Ω\∪N
n=1 Ωn

N
X

∆y G(x, y) − κ20 G(x, y) utm (y)dy +
n=1



Z

G(x, y)
∂Ωn

Dn ∂utm
D0 ∂ν

−

(y) −

∂
G(x, y)utm
∂ν(y)

−


(y) dσ(y).

Using Green’s identity on the last term and by Eq. (16), one arrives at
Z
Z


uim (x) = −
∆y G(x, y) − κ20 G(x, y) utm (y)dy +
G(x, y)∆utm (y) − ∆y G(x, y)utm (y) dy
Ω\∪N
n=1 Ωn

+

N
X

1
D
0
n=1

∪N
n=1 Ωn

Z
G(x, y)δD
∂Ωn

∂utm
∂ν

−

(y)dσ(y).

After fairly easy manipulations, it can be seen that
Z
Z


1
uim (x) = −
∆y G(x, y) − κ20 G(x, y) utm (y)dy +
G(x, y) δµ(y)utm (y) − ∇ · δD(y)∇utm (y) dy
D0 ∪N
Ω
n=1 Ωn
Z
N
t
X
∂u
1
[21]
+
G(x, y)δD m (y)dσ(y).
D
∂ν −
0 ∂Ωn
n=1

n=1

where

RA
FT

Note that the first term on the RHS of Eq. (21) is nothing other than utm thanks to Eq. (17). Moreover, one last
application of the Green’s theorem on the second term furnishes
Z

1
s
um (x) = −
δµ(y)G(x, y)utm (y) + δD(y)∇G(x, y) · ∇utm (y) dy
D0 ∪N
n=1 Ωn
 
Z
δµ
Λm (x, y) ·
dy,
∀x ∈ Ω,
[22]
=
δD
N
∪
Ωn

Λm (x, y) := −D0 −1 G(x, y)utm (y)

>
∇G(x, y) · ∇utm (y) .

D

Eq. (22) is the Lippmann-Schwinger integral equation that relates the scattered field to the perturbations in the
optical properties. In particular, when the absorption perturbations are only present, the scattered photon flux
measurement at the boundary is given by
usm (x) =Mm [f ](x),

where

Mm [f ](x) := −

1
D0

Z
∪N
n=1 Ωn

∀x ∈ ∂Ω,

G(x, y)utm (y)f (y)dy,

 
f := δµ .

B. Bulk optical property estimation
The uniform bulk optical properties were found by fitting the experimental data to the model based data (diffusion
equation in this case) using the iterative Newton-Raphson scheme (46). Under this scheme, two parameters d(log(ρϕ))
dρ
and

d(θ)
dρ

were estimated through linear regression, where ϕ is amplitude, θ is phase of the data (Fig. 9).
O(µa , µ0s ) =



d(log(ρϕmeasured )) d(log(ρϕcalculated ))
−
dρ
dρ

2


+

d(θmeasured ) d(θcalculated )
−
dρ
dρ

2
[23]

The minimization problem in Eq. (23) is the root finding problem which can be efficiently solved by Newton-Raphson
measured
calculated
measured
calculated
))
))
)
)
method. Let α0 = d(log(ρϕdρ
, αc = d(log(ρϕ dρ
, β0 = d(θ dρ
, andβc = d(θ dρ
. Newton-Raphson
0
method find the roots µ = [µa , µs ] of a real-valued function, such that
fα (µ) = αc (µ) − α0
fβ (µ) = βc (µ) − β0
14

|

[24]
Yoo et al.

Fig. 9. Data preprocessing for bulk optical property estimation. The measurements distanced over 51 mm from source to detector were selected, and the bulk optical
properties are estimated using this data by Newton-Raphson scheme.

The method alternates between two functions of Eq. (24) and updates the optical properties using Eq. (25) at each
iterative step.
fα (µi )
fα0 (µi )
fβ (µi )
= µi − 0 i
fβ (µ )

µi=1 = µi −

[25]

RA
FT

µi=1

D

where the derivative f 0 (µi ) is computed numerically.

Yoo et al.

|

December 5, 2017

|

vol. XXX

|

no. XX

|

15

