The decay behavior of specific intensity is studied for spatial-frequency domain imaging (SFDI). It is shown using the radiative transport equation that the decay is given by a superposition of different decay modes, and the decay rates of these modes are determined by spatial frequencies and Case’s eigenvalues. This explains why SFDI can focus on shallow regions. The fact that light with nonzero spatial frequency rapidly decays makes it possible to exclusively extract optical properties of the top layer of a layered medium. We determine optical properties of the top layer of a solid phantom. This measurement is verified with different layered media of numerical phantoms.
1. INTRODUCTION
In near-infrared spectroscopy, light illumination in the spatial-frequency domain has been developed as a tool that is as concise as continuous-wave illumination and as informative as frequency-domain illumination [1]. Spatial-frequency domain imaging (SFDI) is capable of determining both absorption and scattering coefficients from time-independent measurements. SFDI is used mainly to extract optical properties at depths on the order of millimeters. SFDI was used for imaging skin flap oxygenation during reconstructive breast surgery [2]. It was also used to record biochemical compositional changes in a port wine stain after laser therapy [3]. Tissue optical properties of a human volar forearm were estimated by SFDI [4]. Burn wounds were examined by SFDI coupled with laser speckle imaging [5]. See a recent review by Angelo et al. [6] and references therein.
The abovementioned works show that SFDI can exclusively study shallow regions near the skin. In SFDI, spatially modulated incident beams rapidly decay in biological tissue [7]. In this paper, we further investigate this feature of SFDI. This property of SFDI is advantageous when we are interested in measuring optical properties of shallow regions. Even when the thickness is thin, we can assume the half space, which is unbounded in the depth direction. In this paper, we identify optical properties of the top layer of a layered medium. It is not possible to extract optical properties of the top layer in the standard setting of near-infrared spectroscopy, in which optical fibers are attached on the top of a layered medium, because near-infrared light propagates not only in the top layer but reaches deeper layers.
In this paper, the ability to extract optical properties of the top layer is tested by different numerical phantoms. Moreover, optical properties of the top layer of a solid phantom are determined. As a numerical tool for this parameter identification, we demonstrate that the numerical scheme for the radiative transport equation (RTE) based on the method of rotated reference frames [8,9] provides an efficient numerical algorithm for SFDI partially because the method relies on the Fourier transform in the spatial-frequency domain.
Noting that the solution to the RTE is expressed as a superposition of three-dimensional singular eigenfunctions [10], we investigate the asymptotic behavior of the solution. The effect of the spatial frequency ${q_0}$ on the decay of the specific intensity is found. The longest-lived mode is controlled by ${q_0}$ and the largest Case’s eigenvalue. This finding gives a theoretical reason that shallow regions can be exclusively studied by SFDI.
Then we numerically solve the RTE in the half space by the method of rotated reference frames [8]. The method of rotated reference frames in the half space was developed for a point source [9,11] and for a spatially oscillating source [12]. In [9,13], the subtraction of the ballistic term was considered. Using the method of rotated reference frames, the effect of surface scattering in SFDI was studied [14]. The inverse problem is solved by the Levenberg–Marquardt algorithm. With our approach, optical parameters of the top layer of a layered solid phantom made of epoxy resin was determined within 1 s on a laptop computer.
The remainder of this paper is organized as follows. In Section 2, we study the asymptotic behavior of the specific intensity to consider the penetration depth of near-infrared light illuminated by a spatially modulated source. In Section 3, reconstruction of optical properties of the top layer of layered media is considered for different numerical phantoms. Optical properties are estimated using a solid phantom in Section 4. Discussion and conclusions are given in Section 5. In Appendix A, two types of the diffusion approximation are introduced. Appendix B is devoted to the eigenmode expansion and numerical algorithm of the RTE.
2. DECAY BEHAVIOR
Let us consider near-infrared light propagation in the half space. Let $\Omega$ be the half space, i.e.,
(1)$$\begin{split}\Omega = \left\{{\textbf{r} \in {{\mathbb{R}}^3}; - \infty \lt x \lt \infty , - \infty \lt y \lt \infty , 0 \lt z \lt \infty} \right\}.\end{split}$$
Let $\partial \Omega$ be the boundary of $\Omega$, i.e., the $x - y$ plane. The specific intensity at position $\textbf{r}$ ($\textbf{r} = ({\boldsymbol \rho} ,z)^T$, ${\boldsymbol \rho} = (x,y)^T$) in direction $\hat{\textbf{s}}$ is denoted by $I(\textbf{r},\hat{\textbf{s}})$, where $\hat{\textbf{s}}$ is a unit vector specified by the polar angle $\vartheta \in [0,\pi]$ and azimuthal angle $\varphi \in [0,2\pi)$. Let $d\hat{\textbf{s}}$ denote $\sin \vartheta d\vartheta d\varphi$. The RTE is written as (2)$$\left\{{\begin{array}{l}({\hat{\textbf{s}} \cdot \nabla + {\mu _t}})I(\textbf{r},\hat{\textbf{s}}) = {\mu _s}\int_{{{\mathbb S}^2}} p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime)I(\textbf{r},\hat{\textbf{s}}^\prime) {\rm{d}}\hat{\textbf{s}}^\prime ,\\(\textbf{r},\hat{\textbf{s}}) \in \Omega \times {{\mathbb S}^2},\\I(\textbf{r},\hat{\textbf{s}}) = {R_{\mathfrak{n}}}(\hat{\textbf{s}} \cdot \hat{\textbf{z}})I(\textbf{r},{{\hat{\textbf{s}}}_R}) + {I_{{\rm{inc}}}}(\textbf{r},\hat{\textbf{s}}),\quad (\textbf{r},\hat{\textbf{s}}) \in {\Gamma _ -},\end{array}} \right.$$
where total attenuation ${\mu _t}$ is the sum of absorption coefficient ${\mu _a}$ and scattering coefficient ${\mu _s}$, which are both assumed to be positive constants, and $p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime)$ is the scattering phase function. We assume that $p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime)$ is given by (3)$$p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime) = \sum\limits_{l = 0}^{{l_{{\max}}}} \sum\limits_{m = - l}^l {{\rm{g}}^l}{Y_{{lm}}}(\hat{\textbf{s}})Y_{{lm}}^*(\hat{\textbf{s}}^\prime),$$
where ${l_{{\max}}}$ is a positive integer, ${\rm{g}} \in (- 1,1)$ is a constant, and the superscript $*$ denotes a complex conjugate. Spherical harmonics ${Y_{{lm}}}(\hat{\textbf{s}})$ are defined by (4)$${Y_{{lm}}}(\hat{\textbf{s}}) = \sqrt {\frac{{2l + 1}}{{4\pi}}\frac{{(l - m)!}}{{(l + m)!}}} P_l^m(\cos \vartheta){e^{{im\varphi}}},$$
where $P_l^m(\mu)$ are associated Legendre polynomials. Throughout the paper, we set ${l_{{\max}}} = 9$. Moreover, (5)$$\begin{split}{{\Gamma _ -}}&= \left\{{(\textbf{r},\hat{\textbf{s}}) \in \partial \Omega \times {{\mathbb S}^2};\;\nu (\textbf{r}) \cdot \hat{\textbf{s}} \lt 0} \right\}\\&= \left\{{(\textbf{r},\hat{\textbf{s}}) \in \partial \Omega \times {\mathbb S}_ + ^2} \right\},\end{split}$$
where $\nu (\textbf{r})$ is the outer unit vector normal to $\textbf{r} \in \partial \Omega$, and ${\mathbb S}_ + ^2$ denotes the set of unit vectors in inward directions. We give the incident beam ${I_{{\rm{inc}}}}(\textbf{r},\hat{\textbf{s}})$ as (6)$${I_{{\rm{inc}}}}(\textbf{r},\hat{\textbf{s}}) = {e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\delta (\hat{\textbf{s}} - \hat{\textbf{z}}),\quad {\textbf{q}_0} \in {{\mathbb R}^2},$$
where $\hat{\textbf{z}}$ is the unit vector in the positive $z$ direction. The Fresnel reflection for the ratio $\mathfrak{n}$ between the refractive indices inside and outside is also considered in the boundary condition. The direction ${\hat{\textbf{s}}_R}$ is specified by the polar angle $\pi - \vartheta$ and azimuthal angle $\varphi$. Assuming unpolarized light, the Fresnel coefficient ${R_{\mathfrak{n}}}(\mu)$ ($0 \lt \mu\le 1$) is given by [15] (7)$${R_{\mathfrak{n}}}(\mu) = \left\{{\begin{array}{ll}{\frac{1}{2}\left({{{\left({\frac{{\mu - \mathfrak{n}{\mu _0}}}{{\mu + \mathfrak{n}{\mu _0}}}} \right)}^2} + {{\left({\frac{{{\mu _0} - \mathfrak{n}\mu}}{{{\mu _0} + \mathfrak{n}\mu}}} \right)}^2}} \right)}&{\quad {\rm{for}}\;\mu \ge {\mu _c},}\\1&{\quad {\rm{for}}\;\mu \lt {\mu _c},}\end{array}} \right.$$
where ${\mu _0} = \sqrt {1 - {\mathfrak{n}^2}(1 - {\mu ^2})}$, and ${\mu _c} = \sqrt {{\mathfrak{n}^2} - 1} /\mathfrak{n}$.Suppose that ${\textbf{r}_d}$ is a point on $\partial \Omega$. We detect the hemispheric flux
(8)$$\begin{array}{*{20}{l}}{{J_ +}({\textbf{r}_d})}&= \int_0^{2\pi} \int_\pi ^{2\pi} (\cos \vartheta)I({\textbf{r}_d},\hat{\textbf{s}})\sin \vartheta {\rm d}\vartheta {\rm d}\varphi\\&= - {A_{{\rm{RTE}}}}({q_0}){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}},\end{array}$$
where ${q_0} = |{\textbf{q}_0}|$, and ${A_{{\rm{RTE}}}}({q_0})$ is given by (B41) in Appendix B.The fluence of the specific intensity is asymptotically governed by the diffusion equation. Let $u(\textbf{r})$ denote the solution to the diffusion equation. To do the diffusion approximation, we split the specific intensity into two terms (see Appendix A): $I(\textbf{r},\hat{\textbf{s}}) = {I_0}(\textbf{r},\hat{\textbf{s}}) + {I_1}(\textbf{r},\hat{\textbf{s}})$, where ${I_0}(\textbf{r},\hat{\textbf{s}})$ satisfies
(9)$$\left\{{\begin{array}{ll}{\left({\hat{\textbf{s}} \cdot \nabla + \bar\mu} \right){I_0}(\textbf{r},\hat{\textbf{s}}) = 0,}&{\quad (\textbf{r},\hat{\textbf{s}}) \in \Omega \times {{\mathbb S}^2},}\\{{I_0}(\textbf{r},\hat{\textbf{s}}) = {I_{{\rm{inc}}}}(\textbf{r},\hat{\textbf{s}}),}&{\quad (\textbf{r},\hat{\textbf{s}}) \in {\Gamma _ -}.}\end{array}} \right.$$
Here, different choices are possible for $\bar\mu$ [16,17]. Then ${I_1}$ is determined depending on the choice of $\bar\mu$. The ${P_1}$ approximation is made for this ${I_1}$, and eventually we arrive at the diffusion equation. We introduce (10)$${\mu _s^\prime} = (1 - {\rm{g}}){\mu _s},\quad {\mu _*} = {\mu _a} + {\mu _s^\prime},\quad {\mu _{{\rm{eff}}}} = \sqrt {3{\mu _a}{\mu _*}} .$$
Probably the most naive choice is (11)$$\bar\mu= {\mu _a} + {\mu _s}.$$
We call this diffusion approximation DA1. In this case, from (A9), (A11), (A13), and (A16), we obtain (12)$$u(\textbf{r}) = {v_{{\rm{DA1}}}}(z){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}},$$
where (13)$$\begin{array}{*{20}{l}}{{v_{{\rm{DA1}}}}(z)}&= \frac{{3{\mu _s}({\mu _*} + {\rm{g}}{\mu _t})}}{{\mu _t^2 -\mu_{{\rm{eff}}}^2 - q_0^2}}\\&\quad\times \left({\frac{{{\mu _t} + 3{\mu _*}/\zeta}}{{\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + 3{\mu _*}/\zeta}}{e^{- \sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} z}} - {e^{- {\mu _t}z}}} \right).\end{array}$$
Another choice is to set [18] (14)$$\bar\mu= {\mu _*}.$$
We refer to this diffusion approximation as DA2. In this case, from (A9), (A25), and (A26), we have (15)$$u(\textbf{r}) = {v_{{\rm{DA2}}}}(z){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}},$$
where (16)$$\begin{array}{*{20}{l}}{{v_{{\rm{DA2}}}}(z)}&= \frac{{3{{\mu ^\prime_s}}{\mu _*}}}{{\mu _*^2 -\mu_{{\rm{eff}}}^2 - q_0^2}}\\&\quad\times \left({\frac{{{\mu _*} + 3{\mu _*}/\zeta}}{{\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + 3{\mu _*}/\zeta}}{e^{- \sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} z}} - {e^{- {\mu _*}z}}} \right).\end{array}$$
Using the diffusion approximation, the detected light ${J_ +}({\textbf{r}_d})$ can be expressed as (17)$${-}{A_{{\rm{DA1}}}}({q_0}){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\quad {\rm{or}}\quad - {A_{{\rm{DA2}}}}({q_0}){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}},$$
depending on DA1 or DA2. They are given by (A19) or (A29) in Appendix A.In the singular-eigenfunction approach [19], the separation constant $\nu$ is either an eigenvalue ${\nu _j}(M) \gt 1$ ($j = 1, \ldots ,{J^M}$) or in the continuous spectrum $(0,1)$. Although there is only one positive eigenvalue ${\nu _0}$ in the case of isotropic scattering, in general there are multiple eigenvalues. We order them as ${\nu _1} \gt {\nu _2} \gt \cdots \gt {\nu _{{J^M}}} \gt 1$ for each $M$.
Case’s method can be extended to three dimensions [10], and the solution to (2) can be expressed as
(18)$$\begin{split}\!\!\!\!I(\textbf{r},\hat{\textbf{s}}) &={ {e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\sum\limits_{M = - {l_{{\max}}}}^{{l_{{\max}}}}}\\[-5pt]&\quad\times {\left[\sum\limits_{j = 0}^{{J^M} - 1} a_j^M\Psi _{{\nu _j}(M)}^M(\hat{\textbf{s}},{\textbf{q}_0}){e^{- {\mu _t}{{ \hat{k}}_z}({\nu _j}(M){q_0})z/{\nu _j}(M)}}\right.}\\[-5pt]&\quad+{ \left.\int_0^1 {a^M}(\nu)\Psi _\nu ^M(\hat{\textbf{s}},{\textbf{q}_0}){e^{- {\mu _t}{{ \hat{k}}_z}(\nu {q_0})z/\nu}} {\rm d}\nu\vphantom{\sum\limits_{j = 0}^{{J^M} - 1} a_j^M\Psi _{{\nu _j}(M)}^M(\hat{\textbf{s}},{q_0}){e^{- {\mu _t}{{ \hat{k}}_z}({\nu _j}(M){q_0})z/{\nu _j}(M)}}}\right]},\end{split}$$
with coefficients $a_j^M$, ${a^M}(\nu)$. Here, $\Psi _{{\nu _j}(M)}^M(\hat{\textbf{s}},{\textbf{q}_0}),\Psi _\nu ^M(\hat{\textbf{s}},{\textbf{q}_0})$ are three-dimensional singular eigenfunctions introduced in Appendix B.We note that ${\nu _0} = {\nu _1}(0) \gt 1$ is the largest eigenvalue. Let us define
(19)$$\begin{array}{*{20}{l}}{{I_0}(\textbf{r},\hat{\textbf{s}})}&= {I_0}({\boldsymbol \rho} ,z,\hat{\textbf{s}})\\&= {e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}a_0^0\Psi _{{\nu _0}}^0(\hat{\textbf{s}},{\textbf{q}_0}){e^{- {\mu _t}{\hat{k}_z}({\nu _0}{q_0})z/{\nu _0}}}.\end{array}$$
When $z$ is large, the contribution of the mode ${I_0}$ dominates: (20)$$\begin{array}{*{20}{l}}{{{\left\| {I(\cdot ,z, \cdot) - {I_0}(\cdot ,z, \cdot)} \right\|}_{{L^\infty}({{\mathbb R}^2};{L^\infty}({{\mathbb S}^2}))}}}\\= o\left({\exp \left({- z\sqrt {{{\left({\frac{{{\mu _t}}}{{{\nu _0}}}} \right)}^2} + q_0^2}} \right)} \right),\end{array}$$
as $z \to \infty$.In the case of isotropic scattering ($g = 0$) [19], the eigenvalue ${\nu _0}$ satisfies $1 - ({\mu _s}/{\mu _t}){\nu _0}\mathop {\tanh}\nolimits^{- 1} (1/{\nu _0}) = 0$. When ${\mu _a} \ll {\mu _s}$ as is typical in biological tissue, we have $1/{\nu _0} \approx \sqrt {3(1 - {\mu _s}/{\mu _t})}$ and
(21)$$\frac{{{\mu _t}}}{{{\nu _0}}} \approx \sqrt {3{\mu _a}{\mu _t}} .$$
Therefore, we have $I \sim \exp (- z\sqrt {3{\mu _a}{\mu _t} + q_0^2})$.In the general case of $g \ne 0$, we can estimate ${\nu _0}$ using the fact that Case’s eigenvalues are approximately obtained as eigenvalues of a tridiagonal matrix $B(M)$ [see (B29)]. It is found [16] that
(22)$$\frac{{\sqrt {1 + \eta}}}{{\sqrt {3\frac{{{\mu _a}}}{{{\mu _t}}}\left({1 - g\frac{{{\mu _s}}}{{{\mu _t}}}} \right)}}} \le {\nu _0} \le \frac{{1 + \sqrt \eta}}{{\sqrt {3\frac{{{\mu _a}}}{{{\mu _t}}}\left({1 - g\frac{{{\mu _s}}}{{{\mu _t}}}} \right)}}},$$
where (23)$$\eta = \frac{4}{5}\frac{{{\mu _a}}}{{{\mu _a} + {\mu _s}\left({1 - {g^2}} \right)}}.$$
When ${\mu _a}$ is small, we have (24)$$\frac{{{\mu _t}}}{{{\nu _0}}} \approx \sqrt {3{\mu _a}({\mu _t} - g{\mu _s})} \approx \sqrt {3{\mu _a}{{\mu ^\prime_s}}} ,$$
where ${\mu ^\prime _s} = (1 - {\rm{g}}){\mu _s}$. Thus, if ${\mu _a}$ is small, we have in general (25)$$I \sim {e^{- z\sqrt {3{\mu _a}{{\mu ^\prime_s}} + q_0^2}}}.$$
When ${\mu ^\prime _s}$ is large, the asymptotic decay of the diffusion equation is given by $\exp (- z\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2})$ in both diffusion approximations given in (12) and (15). Note that ${\mu _{{\rm{eff}}}} \approx \sqrt {3{\mu _a}{{\mu ^\prime_s}}}$ for small ${\mu _a}$. Since ${\mu _t} \lt {\nu _0} \lt {\mu _{{\rm{eff}}}}$ if $\eta \gt 0$ is taken into account, the specific intensity decays slower than the prediction by the diffusion approximation.If ${q_0}$ is large such that ${\mu _t} \lt \sqrt {\mu _{{\rm{eff}}}^2 + q_0^2}$ or ${\mu _*} \lt \sqrt {\mu _{{\rm{eff}}}^2 + q_0^2}$, the asymptotic decay from the diffusion approximation given in (12) or (15) becomes $\exp (- {\mu _t}z)$ or $\exp (- {\mu _*}z)$, and in either case, the asymptotic behaviors of the RTE and diffusion equation are quite different.
3. SPATIAL-FREQUENCY DOMAIN IMAGING
Although different choices are possible, we modulate the illuminating light in the $x$ direction and give the vector ${\textbf{q}_0}$ as
(26)$${\textbf{q}_0} = (2\pi\! f{,0)^T}.$$
When the sample is illuminated by the source, (27)$${I_{{\rm{inc}}}}(\textbf{r},\hat{\textbf{s}}) = \cos (2\pi\! fx)\delta (\hat{\textbf{s}} - \hat{\textbf{z}}).$$
The measured light is expressed as
(28)$${J_ +}({\textbf{r}_d}) = - A({q_0})\cos (2\pi\! f{x_d}),$$
where ${x_d}$ is the first component of ${\textbf{r}_d}$. Here, $A({q_0})$ depends on ${q_0} = 2\pi\! f$.We use the amplitude $A({q_0})$ to reconstruct optical properties. Let ${N_f}$ be the number of spatial frequencies used for reconstruction. We have
(29)$$f = {f_1}, \ldots ,{f_{{N_f}}}.$$
Correspondingly, we write the forward data as $A(q_0^{(i)})$ ($i = 1, \ldots ,{N_f}$). Let ${\textbf{y}} \in {{\mathbb R}^{{N_f}}}$ be a vector defined as ${\textbf{y}} =$$(A(q_0^{(i)}))$. We can express ${A_{{\rm{RTE}}}}({q_0})$ as ${A_{{\rm{RTE}}}}(q_0^{(i)})$ ($i = 1, \ldots ,{N_f}$). Similarly, ${A_{{\rm{DA1}}}}(q_0^{(i)})$ and ${A_{{\rm{DA2}}}}(q_0^{(i)})$ are introduced. Then computed values are stored in a vector ${\textbf{F}} \in {{\mathbb R}^{{N_f}}}$.Parameters ${\mu _a},{\mu _s^\prime}$ are determined by the Levenberg–Marquardt algorithm [20,21]. To run the inversion algorithm with scaled variables ${\xi _1},{\xi _2}$, we express the optical properties as [22,23]
(30)$${\mu _a} =\mu_a^{(0)}{e^{{\xi _1}}},\quad {\mu ^\prime _s} = {\mu ^\prime _s}^{(0)}{e^{{\xi _2}}},$$
where $\mu_a^{(0)},{\mu ^\prime _s}^{(0)}$ are initial guesses. That is, ${\xi _1} = \ln ({\mu _a}/\mu _a^{(0)})$, ${\xi _2} = \ln ({\mu ^\prime _s}/{\mu ^\prime _s}^{(0)})$.We express ${\textbf{F}} = {\textbf{F}}({\boldsymbol{\xi}})$, where ${\boldsymbol{\xi}} = ({\xi _1},{\xi _2}{)^T}$. We wish to find
(31)$${\boldsymbol{\xi}}_{*} = \mathop{{\arg} \,{\min}}\limits_{\boldsymbol{\xi}} \left| \textbf{y} - \textbf{F}({\boldsymbol{\xi}}) \right|^{2}.$$
Let ${{\boldsymbol{\xi}} ^k}$ be the estimated solution at the $k$th iteration. Starting with the initial guess ${\boldsymbol{\xi}} ^{0} = \textbf{0}$, the solution is given by $\boldsymbol{\xi} _{*} = \mathop{\lim}\nolimits_{k \to \infty} \,\boldsymbol{\xi} ^{k}$. The FORTRAN library MINPACK [24] was used for the numerical calculation of the Levenberg–Marquardt method. Reconstructed values are obtained as (32)$${\mu _a} =\mu_a^{(0)}\exp ({\xi _{*,1}}),\quad {\mu ^\prime _s} = {\mu ^\prime _s}(0)\exp ({\xi _{*,2}}).$$
Below, we will perform several parameter identifications using numerical phantoms.A. Thin Slabs
We perform parameter identification for numerical slab phantoms of size $90 \;{\rm{mm}} \times 90 \;{\rm{mm}} \times L$ as shown in Fig. 1. The thickness $L$ changes from 1 mm to 10 mm. The optical parameters of slabs are set to ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$, ${\mu _s} = 10 \;{\rm{m}}{{\rm{m}}^{- 1}}$, and ${\rm{g}} = 0.9$. Moreover, $\mathfrak{n} = 1$ (vacuum boundary condition). We suppose ${\mu _a},{\mu ^\prime _s}$ are unknown. We set ${N_f} = 2$, and ${f_1} = 0.1$, ${f_2} = 0.2$ (${\rm{m}}{{\rm{m}}^{- 1}}$).
As the forward data, the hemispheric flux is computed by Monte Carlo simulations and stored in ${\textbf{y}}$. In each run, ${10^8}$ photons are launched. The vector ${\textbf{F}} = ({A_{{\rm{RTE}}}}(q_0^{(i)}))$ ($i = 1,2$) is computed from the RTE [see (B41)]. For the inverse problem of parameter identification, initial values are set to $(\mu _a^{(0)},\mu _s^{(0)}) = (0.01 \;{\rm{m}}{{\rm{m}}^{- 1}},10 \;{\rm{m}}{{\rm{m}}^{- 1}})$.
Figure 2 shows estimated ${\mu _a}$ and ${\mu ^\prime _s}$ for the slabs. Estimated values are $({\mu _a} {\rm{m}}{{\rm{m}}^{- 1}},{\mu ^\prime _s} {\rm{m}}{{\rm{m}}^{- 1}}) = (0.18,1.2)$, (0.055, 1.1), (0.033, 1.1), (0.023, 1.0), (0.023, 1.0), (0.024, 1.1), (0.023, 1.0), (0.021, 1.0), (0.023, 1.0), and (0.024, 1.0) for $L = 1 \;{\rm{mm}}$, 2 mm, 3 mm, 4 mm, 5 mm, 6 mm, 7 mm, 8 mm, 9 mm, and 10 mm, respectively. We see that for slabs of thickness larger than 4 mm, optical properties are correctly obtained.
The behavior in Fig. 2 is implied in the decay in (25) derived in Section 2. In the present situation, we have
(33)$$I \sim {e^{- 2\pi fz}}.$$
For $f = 0.1 \;{\rm{m}}{{\rm{m}}^{- 1}}$, we have ${e^{- 2\pi f \cdot 3}} = 0.15$ and ${e^{- 2\pi f \cdot 4}} = 0.08$. Hence, the specific intensity is reduced by more than one tenth when the thickness of the slab is 4 mm or larger. Recently, an intensive study of Monte Carlo look-up tables was reported for relations between the penetration depths of photons and spatial frequencies [25]. Our conclusion in Fig. 2 is consistent with their results.B. Two-Layer Media
Let us consider two-layer media, which have top and bottom layers. As shown in Fig. 3, the thickness of the top layer is 6 mm, and the bottom layer has a thickness of 30 mm, which can be regarded as a semi-infinite medium. We set ${N_f} = 2$, and ${f_1} = 0.1$, ${f_2} = 0.2$ (${\rm{m}}{{\rm{m}}^{- 1}}$). The scattering coefficient and anisotropic factor are fixed to ${\mu _s} = 10 \;{\rm{m}}{{\rm{m}}^{- 1}}$ and ${\rm{g}} = 0.9$, respectively, in the entire medium. The refractive index is set to $\mathfrak{n} = 1.4$. The absorption coefficient in the top layer is ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$. The absorption coefficient ${\mu _a}$ in the bottom layer takes values $0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ and $0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$. We use ${A_{{\rm{RTE}}}}(q_0^{(i)})$ ($i = 1,2$) for reconstruction. The forward data were computed by Monte Carlo simulations. Initial values were set to $(\mu _a^{(0)},\mu _s^{(0)}) = (0.01\;{\rm{m}}{{\rm{m}}^{- 1}},10\;{\rm{m}}{{\rm{m}}^{- 1}})$.
Table 1. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) of the Top Layer of the Two-Layer Medium by RTE for ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Below, we present obtained ${\mu _a},{\mu ^\prime _s}$ in the top layer when the absorption coefficient of the bottom layer is $0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ (Table 1) and $0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$ (Table 2). To see how the estimated value is close to the true value, in Tables 1 and 2, we also give the relative error, which is defined as $|({\rm{estimated \;value}} - {\rm{true \;value}})/{\rm{true \;value}}|$. Numerical results show the reconstructed values are not affected by optical properties of the bottom layer.
Table 2. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) of the Top Layer of the Two-Layer Medium by RTE for ${\mu _a} = 0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
C. Three-Layer Media
Here, we consider more complex media that have more than two layers [26,27]. We reconstruct the optical properties of the top layer of a three-layer medium, shown in Fig. 4. Their optical properties are summarized in Table 3. The depths of layers are, from the top, 10 mm, 2 mm, and 28 mm. From the top, ${\mu _s} = 18$, 3, and 21 (${\rm{m}}{{\rm{m}}^{- 1}}$). In the second layer, ${\mu _a} = 0.004 \;{\rm{m}}{{\rm{m}}^{- 1}}$. The absorption coefficients of the first and third layers vary from $0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ to $0.03\; {\rm{m}}{{\rm{m}}^{- 1}}$. In addition, ${\rm{g}} = 0.9$ and $\mathfrak{n} = 1.4$. We set ${N_f} = 2$, and ${f_1} = 1/15$, ${f_2} = 1/10$ (${\rm{m}}{{\rm{m}}^{- 1}}$). Starting the Levenberg–Marquardt algorithm with initial values $(\mu _a^{(0)},\mu _s^{(0)}) = (0.01\;{\rm{m}}{{\rm{m}}^{- 1}},10\;{\rm{m}}{{\rm{m}}^{- 1}})$, we obtain ${\mu _a},{\mu ^\prime _s}$ after about 10 iterations.
Table 3. Three-Layer Modela
For comparison, reconstructions by the diffusion approximation are also obtained. That is, ${A_{{\rm{DA1}}}}$ in (A19) and ${A_{{\rm{DA2}}}}$ in (A29) are used in addition to ${A_{{\rm{RTE}}}}$ in (8). The forward data were prepared by Monte Carlo simulations.
Tables 4–12 show reconstructed ${\mu _a},{\mu ^\prime _s}$ when the true values in the top layer are ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$, $0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$, and $0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$. In each table, relative errors are also shown. The true value of ${\mu ^\prime _s}$ in the top layer is fixed to $1.8 \;{\rm{m}}{{\rm{m}}^{- 1}}$, as shown in Table 3. We see that reconstructed values by DA1 are closer to those by RTE compared with reconstructed values by DA2. The absorption coefficient of the bottom layer is ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in Tables 4, 7, 10, ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in Tables 5, 8, 11, and ${\mu _a} = 0.03\; {\rm{m}}{{\rm{m}}^{- 1}}$ in Tables 6, 9, 12. In all cases, reconstructed values are not affected by the third layer.
Table 4. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 5. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 6. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.03\; {\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 7. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 8. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.02\; {\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 9. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 10. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 11. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.02 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
Table 12. Reconstructed ${\mu _a}$ and ${\mu ^\prime _s}$ (${\rm{m}}{{\rm{m}}^{- 1}}$) by RTE, DA1, and DA2 When ${\mu _a} = 0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Top Layer and ${\mu _a} = 0.03 \;{\rm{m}}{{\rm{m}}^{- 1}}$ in the Bottom Layer
4. SOLID PHANTOM
Figure 5(a) shows a solid phantom made of epoxy resin. The refractive index of the phantom is 1.58. The phantom has a four-layer structure, and the width of the top layer is about 4 mm. When the phantom was made, optical properties of the top layer were aimed at ${\mu _a} = 0.0165 \;{\rm{m}}{{\rm{m}}^{- 1}}$ and ${\mu ^\prime _s} = 1.3\; {\rm{m}}{{\rm{m}}^{- 1}}$. The setup of the measurement system is shown in Figs. 5(b) and 5(c). As a near-infrared light source, a broadband halogen fiber optic illuminator (Thorlabs, OSL2) with an enhanced infrared replacement bulb (Thorlabs, OSL2BIR) was used. The projected patterns were displayed on a digital micro-mirror device (DMD) module (Keynote Photonics, LC4500 NIR controller). The image was formed on the sample plane through a visible-near-infrared lens (Schneider, large format F-mount lens, focal length of 28 mm, F/2.8). To capture reflection images, a camera with enhanced near-infrared sensitivity (Ximea, MQ-13RG-E2, 1280x1024 pixels, monochrome) with a visible-near-infrared lens (Edmund Optics, C series VIS-NIR fixed focal length lens, focal length of 16 mm, F/1.6) was used. The polarizer and analyzer placed in the crossed nicols configuration were inserted into the measurement system to remove the specular reflection from the sample. A band-pass filter (Edmund Optics, hard coated OD 4.0 25 nm band-pass filter, center wavelength of 800 nm, FWHM of 25 nm) was placed in front of the camera to extract the 800 nm wavelength. The top layer of the phantom was illuminated by the spatially modulated light, and the reflected light was detected by the camera.
The source term in experiments is given by
(34)$${I_{{\rm{inc}}}}(\textbf{r},\hat{\textbf{s}}) = \frac{{{S_0}}}{2}\left[{1 + \cos \!\left({2\pi\! fx + \alpha} \right)} \right]\delta (\hat{\textbf{s}} - \hat{\textbf{z}}),$$
where $\alpha = 2\pi\! p/3$ ($p = 0,1,2$). We set ${N_f} = 2$, and ${f_1} = 0.1$, ${f_2} = 0.2$ (${\rm{m}}{{\rm{m}}^{- 1}}$). We have (35)$${J_ +}({\textbf{r}_d};\alpha) = \frac{{{S_0}}}{2}\left({{M_{{\rm{DC}}}} + {M_{{\rm{AC}}}}\cos (2\pi\! f{x_d} + \alpha)} \right),$$
where ${M_{{\rm{DC}}}},{M_{{\rm{AC}}}}$ depend on ${\textbf{r}_d}$ and $f$ in general. Let us write $J_ + ^{(p)} = {J_ +}({\textbf{r}_d};2\pi\! p/3)$. By a straightforward calculation, we have (36)$$\begin{array}{*{20}{l}}{{{({{M_{{\rm{AC}}}}} )}^2} = \frac{2}{9}}\\\quad\times \big[{{{\big({J_ + ^{(0)} - J_ + ^{(1)}} \big)}^2} + {{\big({J_ + ^{(1)} - J_ + ^{(2)}} \big)}^2} + {{\big({J_ + ^{(2)} - J_ + ^{(0)}} \big)}^2}} \big].\end{array}$$
Moreover, we can write [7] (37)$${M_{{\rm{AC}}}} = {M^{{\exp}}}({\textbf{r}_d},f){A^{{\exp}}}(q_0^{(i)}),$$
where ${M^{{\exp}}}$ (${\gt}0$) is a constant determined by the optical system, and $q_0^{(i)} = {f_i}$ ($i = 1,2$).Let us consider ${M_{{\rm{AC,}} {\rm{ref}}}}$ of a reference medium, for which $A_{{\rm{ref}}}^{{\exp}}({q_0})$ can be numerically computed. Then we have [7]
(38)$${A^{{\exp}}}({q_0}) = \left\langle {\frac{{{M_{{\rm{AC}}}}}}{{{M_{{\rm{AC,}} {\rm{ref}}}}}}} \right\rangle A_{{\rm{ref}}}^{{\exp}}({q_0}),$$
where $\langle \cdot \rangle$ means the average in space.We use the bottom layer of the solid phantom as the reference medium. Since the fourth layer of the phantom has the thickness 47 mm, it can be regarded as the half space even for time-resolved measurements in which two optical fibers are vertically attached on the bottom side of the phantom. From time-resolved measurements by TRS-80 (Hamamatsu Photonics), we found ${\mu _a} = 0.010 \;{\rm{m}}{{\rm{m}}^{- 1}}$ and ${\mu^\prime_s} = 1.4 \;{\rm{m}}{{\rm{m}}^{- 1}}$ for the bottom layer. With these optical parameters, $A_{{\rm{ref}}}^{{\exp}}(q_0^{(i)})$ can be computed from the RTE as ${A_{{\rm{RTE}}}}(q_0^{(i)})$. Thus, ${A^{{\exp}}}(q_0^{(i)})$ in (38) are prepared and stored in a vector ${\textbf{y}} \in {{\mathbb R}^{{N_f}}}$.
The Levenberg–Marquardt algorithm was run with initial values $\mu_a^{(0)} = 0.01 \;{\rm{m}}{{\rm{m}}^{- 1}}$ and ${\mu _s^{(0)}} = 10 \;{\rm{m}}{{\rm{m}}^{- 1}}$. We obtain
(39)$${\mu _a} = 0.016 \;{\rm{m}}{{\rm{m}}^{- 1}},\quad {\mu ^\prime _s} = 1.0 \;{\rm{m}}{{\rm{m}}^{- 1}}.$$
Other choices of initial guess, for example, $\mu_a^{(0)} = 0.02\;{\rm{m}}{{\rm{m}}^{- 1}}$ and ${\mu ^\prime _s}^{(0)} = 2\;{\rm{m}}{{\rm{m}}^{- 1}}$, give the same ${\mu _a},{\mu ^\prime _s}$ given in (39). The computation time was less than 1 s with a laptop computer (MacBook Pro with 2.3 GHz Intel Core i5 and 8 GB memory).5. DISCUSSION AND CONCLUSION
Taking advantage of the fact that near-infrared light decays rapidly for nonzero spatial frequencies, in this paper, we estimated optical properties of the top layer of the layered phantom. Indeed, SFDI has been used for the parameter identification of the top layer. Various numerical experiments developed in the present paper confirm that results by SFDI are not affected by deeper layers. Figure 2 indicates that the necessary width of the top layer is 4 mm (when $1/{\mu ^\prime _s} = 1 \;{\rm{mm}}$). Our investigation in the asymptotic limit provides a theoretical reason for the results of numerical experiments. Using the numerical algorithm based on the method of rotated reference frames, optical properties of the top layer of a layered solid phantom were also determined.
Our numerical results show that ${\mu ^\prime _s}$ is more accurately reconstructed than ${\mu _a}$. One of the reasons for this inaccuracy is the ill-posedness of the parameter identification. Future study on regularization might improve the results.
When two optical fibers are attached to biological tissue in the direction perpendicular to its surface with a separation of a few centimeters, the detected reflected near-infrared light contains photons from depths more than a centimeter. Although this very feature makes it possible to study the function of the human brain through the neurovascular coupling [28], the detected light is affected by optical properties of different layers. In particular, the signal from the brain is affected by skin blood flow in the scalp [29]. In this conventional way, it is not possible to extract only information on shallow regions. Photons that travel deep inside biological tissue can be excluded if the separation of two optical fibers is reduced. However, then measurements have to be conducted in a tiny space, and other difficulties related to measurements arise (see [30] and references therein). The SFDI measurement setup described in Fig. 5(a) is free from such difficulties.
For our approach to work, the top layer of a layered random medium has to be regarded as a semi-infinite medium. The necessary width of the top layer depends on spatial frequency. This can be checked by a test parameter identification for numerical phantoms.
When ${\mu _a}$ is not small, the decay of the specific intensity deviates from the diffusive decay and is given by (20) because then ${\mu _t}/{\nu _0}$ cannot be approximated by $\sqrt {3{\mu _a}{{\mu ^\prime_s}}}$. Moreover if the depth is not large compared to ${\nu _0}/{\mu _t}$ or $1/{q_0}$, not only ${I_0}$ but other modes contribute, and the decay is given by the superposition of different decays shown in (18).
APPENDIX A: DIFFUSION APPROXIMATION
We begin by decomposing $I$ into the following two terms:
(A1)$$I(\textbf{r},\hat{\textbf{s}}) = {I_b}(\textbf{r},\hat{\textbf{s}}) + {I_s}(\textbf{r},\hat{\textbf{s}}).$$
The ballistic term
${I_b}$ and scattering term
${I_s}$ satisfy
(A2)$$\left\{{\begin{array}{ll}{\left({\hat{\textbf{s}} \cdot \nabla + {\mu _a} + {\mu _s}} \right){I_b}(\textbf{r},\hat{\textbf{s}}) = 0,}&{\quad (\textbf{r},\hat{\textbf{s}}) \in \Omega \times {{\mathbb S}^2},}\\{{I_b}(\textbf{r},\hat{\textbf{s}}) = {I_{{\rm{inc}}}}(\textbf{r},\hat{\textbf{s}}),}&{\quad (\textbf{r},\hat{\textbf{s}}) \in {\Gamma _ -}}\end{array}} \right.$$
and
(A3)$$\left\{{\begin{array}{ll}{\left({\hat{\textbf{s}} \cdot \nabla + {\mu _a} + {\mu _s}} \right){I_s}(\textbf{r},\hat{\textbf{s}}) = {\mu _s}\int_{{{\mathbb S}^2}} p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime){I_s}(\textbf{r},\hat{\textbf{s}}^\prime) {\rm{d}}\hat{\textbf{s}}^\prime}\\{\quad \quad + S(\textbf{r},\hat{\textbf{s}}),\quad (\textbf{r},\hat{\textbf{s}}) \in \Omega \times {{\mathbb S}^2},}\\{{I_s}(\textbf{r},\hat{\textbf{s}}) = {R_{\mathfrak{n}}}(\hat{\textbf{s}} \cdot \hat{\textbf{z}}){I_s}(\textbf{r},{{\hat{\textbf{s}}}_R}),\quad (\textbf{r},\hat{\textbf{s}}) \in {\Gamma _ -}.}\end{array}} \right.$$
Here, the source term for
${I_s}$ is given by
(A4)$$S(\textbf{r},\hat{\textbf{s}}) = {\mu _s}\int_{{{\mathbb S}^2}} p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime){I_b}(\textbf{r},\hat{\textbf{s}}^\prime) {\rm{d}}\hat{\textbf{s}}^\prime .$$
Since
(A5)$${I_b}(\textbf{r},\hat{\textbf{s}}) = {e^{- {\mu _t}z}}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\delta (\hat{\textbf{s}} - \hat{\textbf{z}}),$$
we have
(A6)$$S(\textbf{r},\hat{\textbf{s}}) = {\mu _s}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}{e^{- {\mu _t}z}}p(\hat{\textbf{s}},\hat{\textbf{z}}).$$
Suppose that
${I_s}$ weakly depends on
$\hat{\textbf{s}}$ and can be written as
(A7)$${I_s}(\textbf{r},\hat{\textbf{s}}) = \frac{1}{{4\pi}}u(\textbf{r}) + \frac{3}{{4\pi}}\textbf{J}(\textbf{r}) \cdot \hat{\textbf{s}},$$
where
(A8)$$u(\textbf{r}) = \int_{{{\mathbb S}^2}} {I_s}(\textbf{r},\hat{\textbf{s}}) {\rm d}\hat{\textbf{s}},\quad \textbf{J}(\textbf{r}) = \int_{{{\mathbb S}^2}} \hat{\textbf{s}}{I_s}(\textbf{r},\hat{\textbf{s}}) {\rm{d}}\hat{\textbf{s}}.$$
Let us write
(A9)$$u(\textbf{r}) = {v_{{\rm{DA1}}}}(z){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}.$$
Then from (
A3), we obtain
(A10)$$\frac{{{\partial ^2}}}{{\partial {z^2}}}{v_{{\rm{DA1}}}} - \left({\mu _{{\rm{eff}}}^2 + q_0^2} \right){v_{{\rm{DA1}}}} = - B{e^{- {\mu _t}z}},$$
where
(A11)$$B = 3{\mu _*}{\mu _s}\left({1 + {\rm{g}}\frac{{{\mu _t}}}{{{\mu _*}}}} \right).$$
We note that
(A12)$$\textbf{J}(\textbf{r}) = - \frac{1}{{3{\mu _*}}}\nabla u(\textbf{r}) + \frac{{g{\mu _s}}}{{{\mu _*}}}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}{e^{- {\mu _t}z}}\hat{\textbf{z}}.$$
We obtain
(A13)$${v_{{\rm{DA1}}}}(z) = \frac{{- B}}{{\mu _t^2 -\mu_{{\rm{eff}}}^2 - q_0^2}}{e^{- {\mu _t}z}} + {C_1}{e^{- \sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} z}},$$
where
${C_1}$ is determined from the boundary condition. If we assume the diffuse boundary condition such that
(A14)$${-}{D_0}\frac{{\partial u}}{{\partial z}} + \frac{1}{\zeta}u = 0,$$
where
${D_0} = 1/(3{\mu _*})$, and
$\zeta = 2(1 + {r_d})/(1 - {r_d})$, with
(A15)$$\!\!\!{r_d} = - 1.4399{\mathfrak{n}^{- 2}} + 0.7099{\mathfrak{n}^{- 1}} + 0.6681 + 0.0636\mathfrak{n},$$
we obtain
(A16)$${C_1} = \frac{B}{{\mu _t^2 -\mu_{{\rm{eff}}}^2 - q_0^2}}\frac{{\zeta {D_0}{\mu _t} + 1}}{{\zeta {D_0}\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + 1}}.$$
On the boundary at
${\textbf{r}_d}$, we have
(A17)$$\!\!\!\!u({\textbf{r}_d}) = \frac{{3{\mu _*}{\mu _s}(1 + {\rm{g}}{\mu _t}/{\mu _*}){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}}}{{\left({\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + {\mu _t}} \right)\left({\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + 3{\mu _*}/\zeta} \right)}}.$$
We obtain
(A18)$$\begin{split}{{J_ +}({\textbf{r}_d})}&= \int_{{\mathbb S}_ - ^2} (\cos \vartheta)\left({\frac{1}{{4\pi}}u({\textbf{r}_d}) + \frac{3}{{4\pi}}\textbf{J}({\textbf{r}_d}) \cdot \hat{\textbf{s}}} \right) {\rm d}\hat{\textbf{s}}\\&= - \frac{1}{4}u({\textbf{r}_d}) + \frac{1}{2}{J_z}({\textbf{r}_d})\\&= - {A_{{\rm{DA1}}}}({{q}_0}){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}},\end{split}$$
where
$u({\textbf{r}_d})$ is given in (
A17). Here, we introduced
(A19)$${A_{{\rm{DA1}}}}({q_0}) = \left({\frac{1}{4} + \frac{1}{{2\zeta}}} \right){v_{{\rm{DA1}}}}(0) - \frac{{{\rm{g}}{\mu _s}}}{{2{\mu _*}}}.$$
Next, instead of (
A2), let us introduce the ballistic term as
$(\hat{\textbf{s}} \cdot \nabla + {\mu _*}){I_0} = 0$; we obtain
(A20)$${I_0}(\textbf{r},\hat{\textbf{s}}) = {e^{- {\mu _*}z}}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\delta (\hat{\textbf{s}} - \hat{\textbf{z}}),$$
and
(A21)$$\begin{split}{S(\textbf{r},\hat{\textbf{s}})}&= ({\mu _*} - {\mu _t}){I_0}(\textbf{r},\hat{\textbf{s}}) + {\mu _s}\int_{{{\mathbb S}^2}} p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime){I_0}(\textbf{r},\hat{\textbf{s}}^\prime) {\rm d}\hat{\textbf{s}}^\prime\\ &= \left({- {\mu _s}{\rm{g}}\delta (\hat{\textbf{s}} - \hat{\textbf{z}}) + {\mu _s}p(\hat{\textbf{s}},\hat{\textbf{z}})} \right){e^{- {\mu _*}z}}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}.\end{split}$$
We can write
(A22)$$u(\textbf{r}) = {v_{{\rm{DA2}}}}(z){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}.$$
After similar calculations, we obtain
(A23)$$\textbf{J}(\textbf{r}) = - {D_0}\nabla u(\textbf{r}),$$
and
(A24)$$\frac{{{\partial ^2}}}{{\partial {z^2}}}{v_{{\rm{DA2}}}} - \left({\mu _{{\rm{eff}}}^2 + q_0^2} \right){v_{{\rm{DA2}}}} = - 3{\mu _*}{\mu ^\prime _s}{e^{- {\mu _*}z}}.$$
Hence,
(A25)$${v_{{\rm{DA2}}}}(z) = \frac{{- 3{\mu _*}{{\mu ^\prime_s}}}}{{\mu _*^2 -\mu_{{\rm{eff}}}^2 - q_0^2}}{e^{- {\mu _*}z}} + {C_2}{e^{- \sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} z}},$$
where
(A26)$${C_2} = \frac{{3{\mu _*}{{\mu ^\prime_s}}}}{{\mu _*^2 -\mu_{{\rm{eff}}}^2 - q_0^2}}\frac{{\zeta {D_0}{\mu _*} + 1}}{{\zeta {D_0}\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + 1}}.$$
We obtain
(A27)$$u({\textbf{r}_d}) = \frac{{3{\mu _*}{{\mu ^\prime_s}}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}}}{{\left({\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + {\mu _*}} \right)\left({\sqrt {\mu _{{\rm{eff}}}^2 + q_0^2} + 3{\mu _*}/\zeta} \right)}}.$$
In this case, we obtain
(A28)$${J_ +}({\textbf{r}_d}) = - {A_{{\rm{DA2}}}}({q_0}){e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}},$$
where
(A29)$${A_{{\rm{DA2}}}}({q_0}) = \left({\frac{1}{4} + \frac{1}{{2\zeta}}} \right){v_{{\rm{DA2}}}}(0).$$
The difference between (
A17) and (
A27) is small when
$g$ is small. In this paper, we use (
A17), for which the decomposition is compatible with the ballistic subtraction developed in Appendix
B.
APPENDIX B: THE METHOD OF ROTATED REFERENCE FRAMES
For the method of rotated reference frames in the half space, the subtraction of the ballistic term was considered [13]. Here, we will compute the hemispheric flux ${J_ +}$ following [13] with a spatially oscillating source term.
1. PRELIMINARY
We begin with the one-dimensional RTE:
(B1)$$\left\{{\begin{array}{rr}{\left({\cos \vartheta \frac{\partial}{{\partial z}} + {\mu _t}} \right){I_1}(z,\hat{\textbf{s}}) = {\mu _s}\int_{{{\mathbb S}^2}} p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime){I_1}(z,\hat{\textbf{s}}^\prime) {\rm{d}}\hat{\textbf{s}}^\prime ,}\\{(z,\hat{\textbf{s}}) \in (0,\infty) \times {{\mathbb S}^2},}\\\quad{{I_1}(z,\hat{\textbf{s}}) = {R_{\mathfrak{n}}}(\mu){I_1}(z,{{\hat{\textbf{s}}}_R}) + {g_1}(z,\hat{\textbf{s}}),}\\{\hat{\textbf{s}} \in {\mathbb S}_ + ^2\quad {\rm{at}}\;z = 0}\end{array}} \right.,$$
with a source term
${g_1}(z,\hat{\textbf{s}})$. In one-dimensional transport theory, it is known that the solution
${I_1}$ is expressed as [
19]
(B2)$$\begin{split}{{I_1}(z,\hat{\textbf{s}})}&= \sum\limits_{M = - {l_{{\max}}}}^{{l_{{\max}}}} \left[\sum\limits_{j = 1}^{{J^M}} \tilde a_j^M\Phi _{{\nu _j}(M)}^M(\hat{\textbf{s}}){e^{- {\mu _t}z/{\nu _j}(M)}}\right.\\&\quad+ \left.\int_0^1 {{\tilde a}^M}(\nu)\Phi _\nu ^M(\hat{\textbf{s}}){e^{- {\mu _t}z/\nu}} {\rm{d}}\nu\vphantom{\sum\limits_{j = 1}^{{J^M}} \tilde a_j^M\Phi _{{\nu _j}(M)}^M(\hat{\textbf{s}}){e^{- {\mu _t}z/{\nu _j}(M)}}}\right],\end{split}$$
where coefficients
$\tilde a_j^M$,
${\tilde a^M}(\nu)$ are determined from the boundary condition. Here,
$\Phi _\nu ^M(\hat{\textbf{s}})$ (
$\nu = {\nu _j}(M)$ or
$\nu \in (0,1)$) is given by [
19,
31–
33]
(B3)$$\Phi _\nu ^M(\hat{\textbf{s}}) = {\phi ^M}(\nu ,\cos \vartheta){\left({1 - \mathop {\cos}\nolimits^2 \vartheta} \right)^{|M|/2}}{e^{{iM\varphi}}},$$
where
${\phi ^M}(\nu ,\cos \vartheta)$ is called Case’s singular eigenfunction.
Let us introduce ${\sigma _l} \gt 0$ as
(B4)$${\sigma _l} = {\mu _t} - {\mu _s}{{\rm{g}}^l} = {\mu _a} + (1 - {{\rm{g}}^l}){\mu _s}$$
and introduce
${b_l}(m)$ as
(B5)$${b_l}(m) = \sqrt {({l^2} - {m^2})/((4{l^2} - 1){\sigma _l}{\sigma _{l - 1}})} .$$
We consider the normalized Chandrasekhar polynomial
$g_l^m(x)$, which satisfies the following three-term recurrence relation [
34,
35]:
(B6)$$\begin{split}\!\!\!{\frac{x}{{{\mu _t}}}\sqrt {(2l + 1){\sigma _l}} g_l^m(x)}&= {b_{l + 1}}(m)\sqrt {(2l + 3){\sigma _{l + 1}}} g_{l + 1}^m(x)\!\!\!\\[-5pt] &\quad+ {b_l}(m)\sqrt {(2l - 1){\sigma _{l - 1}}} g_{l - 1}^m(x)\end{split}$$
for
$l \gt m$ and
$m \ge 0$. We have
(B7)$$\begin{split}&{g_m^m(x) = \frac{{(2m - 1)!!}}{{\sqrt {(2m)!}}} = \frac{{\sqrt {(2m)!}}}{{{2^m}m!}},}\\&{g_l^{- m}(x) = (- {{1)}^m}g_l^m(x),}\\&{g_l^m(- x) = (- {{1)}^{l + m}}g_l^m(x).}\end{split}$$
Now, eigenvalues
${\nu _j}(M)$ are zeros of
$g_l^M$ as
$l \to \infty$ [
36]. In (
B28), we numerically obtain
${\nu _j}(M)$ and
$\nu \in (0,1)$ as eigenvalues of a tridiagonal matrix.
Let $f(\hat{\textbf{s}})$ be a function of $\hat{\textbf{s}} \in {{\mathbb S}^2}$, which can be expressed as $f(\hat{\textbf{s}}) = \sum\nolimits_{l = 0}^\infty \sum\nolimits_{m = - l}^l {f_{{lm}}}{Y_{{lm}}}(\hat{\textbf{s}})$, with coefficients ${f_{{lm}}}$. We introduce $ {{{\cal R}_{\hat{\textbf{k}}}}}$ for a unit vector $\hat{\textbf{k}} \in {\mathbb C}$ ($\hat{\textbf{k}} \cdot \hat{\textbf{k}} = 1$) as
(B8)$$\mathop {{{\cal R}_{\hat{\textbf{k}}}}} f(\hat{\textbf{s}}) = \sum\limits_{l = 0}^\infty \sum\limits_{m = - l}^l {f_{{lm}}}{Y_{{lm}}}(\hat{\textbf{s}};\hat{\textbf{k}}),$$
where [
8]
(B9)$${Y_{{lm}}}(\hat{\textbf{s}};\hat{\textbf{k}}) = \sum\limits_{m^\prime = - l}^l {e^{- im^\prime {\varphi _{\hat{\textbf{k}}}}}}d_{m^\prime m}^l({\vartheta _{\hat{\textbf{k}}}}){Y_{{lm^\prime}}}(\hat{\textbf{s}}).$$
Here,
${\varphi _{\hat{\textbf{k}}}},{\vartheta _{\hat{\textbf{k}}}}$ are the azimuthal and polar angles of
$\hat{\textbf{k}}$, and
$d_{m^\prime m}^l$ are the Wigner
$d$ matrices. We choose the branch cut of the square root function from zero to
$\infty$, so
$0 \le {{\arg}} (\sqrt z) \lt \pi$ for arbitrary
$z \in {\mathbb C}$. That is, by
$ {{{\cal R}_{\hat{\textbf{k}}}}} $, we measure angles in
$f(\hat{\textbf{s}})$ in the reference frame rotated so that the
$z$ axis lies in the direction of
$\hat{\textbf{k}}$.
Indeed, $\nu$ is either an eigenvalue or in the continuous spectrum. We take the specific form of $\hat{\textbf{k}} = \hat{\textbf{k}}(\nu ,\textbf{q})$ given below, which depends on $\nu$ and $\textbf{q} \in {{\mathbb R}^2}$:
(B10)$$\hat{\textbf{k}}(\nu ,\textbf{q}) = \left({- i\nu \frac{\textbf{q}}{{{\mu _t}}},\;{\hat{\textbf{k}}_z}(\nu q)} \right),\quad \hat{{k}_z}(\nu q) = \sqrt {1 + {{(\nu q/{\mu _t})}^2}} ,$$
where
$q = |\textbf{q}|$. We note that
$\hat{\textbf{k}}(\nu ,\textbf{q}) \cdot \hat{\textbf{k}}(\nu ,\textbf{q}) = 1$. We obtain
(B11)$${\varphi _{\hat{\textbf{k}}}} = \left\{{\begin{array}{ll}{{\varphi _\textbf{q}} + \pi}&{\quad {\rm{for}}\;\nu \gt 0,}\\{{\varphi _\textbf{q}}}&{\quad {\rm{for}}\;\nu \lt 0,}\end{array}} \right.$$
where
${\varphi _\textbf{q}}$ is the polar angle of
$\textbf{q}$, and
(B12)$$\!\!\!\cos {\vartheta _{\hat{\textbf{k}}}} = \hat{\textbf{k}} \cdot \hat{\textbf{z}} = \hat{k_z},\;\; \sin {\vartheta _{\hat{\textbf{k}}}} = \sqrt {1 - \mathop {\cos}\nolimits^2 {\vartheta _{\hat{\textbf{k}}}}} = i|\nu q|.\!$$
We note that
$d_{m^\prime m}^l({\vartheta _{\hat{\textbf{k}}}})$ depends on
$q$ but is independent of
$\varphi _{\textbf{q}}$. Hence, we write
(B13)$$d_{m^\prime m}^l({\vartheta _{\hat{\textbf{k}}}}) = d_{m^\prime m}^l[i\tau (\nu q)].$$
We define
(B14)$$\begin{split}&{\Psi _{{\nu _j}(M)}^M(\hat{\textbf{s}},{\textbf{q}_0})} = {{{\cal R}_{\hat{\textbf{k}}({\nu _j}(M),{\textbf{q}_0})}}} \Phi _{{\nu _j}(M)}^M(\hat{\textbf{s}}),\\[-5pt]&\Psi _\nu ^M(\hat{\textbf{s}},{\textbf{q}_0}) = {{{\cal R}_{\hat{\textbf{k}}(\nu ,{\textbf{q}_0})}}} \Phi _\nu ^M(\hat{\textbf{s}}).\end{split}$$
2. BALLISTIC SUBTRACTION
As was done in Appendix A, we consider ${I_s}$ by subtracting ${I_b}$ from $I$. Let us introduce the particular solution ${I_p}$ as
(B15)$$\begin{split} \left({\hat{\textbf{s}} \cdot \nabla + {\mu _a} + {\mu _s}} \right){I_p}(\textbf{r},\hat{\textbf{s}})&= {\mu _s}\int_{{{\mathbb S}^2}} p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime){I_p}(\textbf{r},\hat{\textbf{s}}^\prime) {\rm{d}}\hat{\textbf{s}}^\prime\\[-2pt] &\quad+ \Theta (z)S(\textbf{r},\hat{\textbf{s}}),\; (\textbf{r},\hat{\textbf{s}}) \in {{\mathbb R}^3} \times {{\mathbb S}^2},\end{split}$$
where
$\Theta (\cdot)$ is the step function. Then we can calculate
${I_s}$ as
(B16)$${I_s} = {I_p} + \psi ,$$
where
$\psi$ satisfies
(B17)$$\left\{{\begin{array}{r}{\left({\hat{\textbf{s}} \cdot \nabla + {\mu _a} + {\mu _s}} \right)\psi (\textbf{r},\hat{\textbf{s}}) = {\mu _s}\int_{{{\mathbb S}^2}} p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime)\psi (\textbf{r},\hat{\textbf{s}}^\prime) {\rm{d}}\hat{\textbf{s}}^\prime ,}\\{(\textbf{r},\hat{\textbf{s}}) \in \Omega \times {{\mathbb S}^2},}\\{\psi (\textbf{r},\hat{\textbf{s}}) = {R_{\mathfrak{n}}}(\hat{\textbf{s}} \cdot \hat{\textbf{z}})\psi (\textbf{r},{{\hat{\textbf{s}}}_R}) + I_{{\rm{inc}}}^{(p)}(\textbf{r},\hat{\textbf{s}}),}\\{(\textbf{r},\hat{\textbf{s}}) \in {\Gamma _ -}}\end{array}} \right.$$
with
(B18)$$I_{{\rm{inc}}}^{(p)}(\textbf{r},\hat{\textbf{s}}) = {R_{\mathfrak{n}}}(\hat{\textbf{s}} \cdot \hat{\textbf{z}}){I_p}(\textbf{r},{\hat{\textbf{s}}_R}) - {I_p}(\textbf{r},\hat{\textbf{s}}),\quad r \in \partial \Omega .$$
3. PARTICULAR SOLUTION
To find ${I_p}(\textbf{r},\hat{\textbf{s}})$, we write
(B19)$${I_p}(\textbf{r},\hat{\textbf{s}}) = {\mu _s}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}{e^{- {\mu _t}z}}\Theta (z)\sum\limits_{l = 0}^{{l_{{\max}}}} \sum\limits_{m = - l}^l {\eta _{{lm}}}{e^{- im{\varphi _{{\textbf{q}_0}}}}}{Y_{{lm}}}(\hat{\textbf{s}}).$$
By multiplying
$Y_{{lm}}^*(\hat{\textbf{s}})$ on both sides of the RTE for
${I_p}$ and integrating over
$\hat{\textbf{s}} \in {{\mathbb S}^2}$, we arrive at the following linear system that determines
${\eta _{{lm}}}$:
(B20)$$\begin{split}&\!\!\!\!{\sum\limits_{l^\prime m^\prime} \left[\frac{{i{q_0}}}{2}\left(- {\delta _{m^\prime ,m - 1}}{\delta _{l^\prime ,l - 1}}\sqrt {(l + m)(l + m - 1)}\right.\right.}\\[-5pt]&\quad+ {\delta _{m^\prime ,m - 1}}{\delta _{l^\prime ,l + 1}}\sqrt {(l^\prime - m^\prime)(l^\prime - m^\prime - 1)}\\[-5pt]&\quad+ {\delta _{m^\prime ,m + 1}}{\delta _{l^\prime ,l - 1}}\sqrt {(l - m)(l - m - 1)}\\[-5pt]&\quad- \left.{\delta _{m^\prime ,m + 1}}{\delta _{l^\prime ,l + 1}}\sqrt {(l^\prime + m^\prime)(l^\prime + m^\prime - 1)}\right)\\[-5pt]&\quad- {\mu _t}{\delta _{m^\prime m}}\left({{\delta _{l^\prime ,l - 1}}\sqrt {{l^2} - {m^2}} + {\delta _{l^\prime ,l + 1}}\sqrt {{{(l^\prime)}^2} - {m^2}}} \right)\\[-5pt]&\quad+\left. {\delta _{m^\prime m}}{\delta _{l^\prime l}}(2l + 1){\sigma _l}\vphantom{\frac{{i{q_0}}}{2}}\right]\frac{{{\eta _{l^\prime m^\prime}}}}{{\sqrt {2l^\prime + 1}}}\\[-5pt]&\quad= {\delta _{m0}}\frac{{(2l + 1){g^l}}}{{\sqrt {4\pi}}}.\end{split}$$
Suppose that the light in direction
$\hat{\textbf{s}} \in {\mathbb S}_ - ^2$ is detected at
${\textbf{r}_d} \in \partial \Omega$. Here,
${\mathbb S}_ - ^2$ denotes the set of unit vectors in outgoing directions. We have
(B21)$${I_p}({\textbf{r}_d},\hat{\textbf{s}}) = {\mu _s}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\sum\limits_{l = 0}^{{l_{{\max}}}} \sum\limits_{m = - l}^l {\eta _{{lm}}}{e^{- im{\varphi _{{\textbf{q}_0}}}}}{Y_{{lm}}}(\hat{\textbf{s}}).$$
4. GENERAL SOLUTION
Since the scattering phase function $p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime)$ depends only on $\hat{\textbf{s}} \cdot \hat{\textbf{s}}^\prime $, we can rewrite (3) as
(B22)$$p(\hat{\textbf{s}},\hat{\textbf{s}}^\prime) = \sum\limits_{l = 0}^{{l_{{\max}}}} \sum\limits_{m = - l}^l {{\rm{g}}^l}{Y_{{lm}}}(\hat{\textbf{s}};\hat{\textbf{k}})Y_{{lm}}^*(\hat{\textbf{s}}^\prime ;\hat{\textbf{k}})$$
for arbitrary
$\hat{\textbf{k}} = \hat{\textbf{k}}(\nu ,\textbf{q})$. We note that
(B23)$$\begin{split} {Y_{{lM}}^*(\hat{\textbf{s}};\hat{\textbf{k}})}&= {{{\cal R}_{\hat{\textbf{k}}}}} Y_{{lM}}^*(\hat{\textbf{s}})\\ &= \sum\limits_{m = - l}^l {e^{im{\varphi _{\hat{\textbf{k}}}}}}d_{{mM}}^l({\vartheta _{\hat{\textbf{k}}}})Y_{{lm}}^*(\hat{\textbf{s}}).\end{split}$$
Let us express the eigenmodes as
(B24)$$\begin{split}&{{\psi _\nu}(\textbf{r},\hat{\textbf{s}},\textbf{q})}\\&\quad= \sum\limits_{l = 0}^{{l_{{\max}}}} \frac{1}{{\sqrt {{\sigma _l}}}}\left\langle {l|{\phi _n}(M)} \right\rangle {Y_{{lM}}}(\hat{\textbf{s}};\hat{\textbf{k}}){e^{- {\mu _t}\hat{\textbf{k}} \cdot \textbf{r}/\nu}}\\&\quad= {e^{i\textbf{q} \cdot {\boldsymbol \rho}}}{e^{- {\mu _t}{\hat{\textbf{k}}_z}(\nu q)z/\nu}}\\&\qquad\times \sum\limits_{l = 0}^{{l_{{\max}}}} \frac{{\left\langle {l|{\phi _\nu}} \right\rangle}}{{\sqrt {{\sigma _l}}}}\sum\limits_{m = - l}^l {e^{- im{\varphi _{\hat{\textbf{k}}}}}}d_{{mM}}^l({\vartheta _{\hat{\textbf{k}}}}){Y_{{lm}}}(\hat{\textbf{s}}).\end{split}$$
We substitute the above
${\psi _\nu}(\textbf{r},\hat{\textbf{s}},\textbf{q})$ in the homogeneous equation of the RTE. By using
$l^\prime = 0, \ldots ,{l_{{\max}}}$ and
$m^\prime = 0, \pm 1, \ldots , \pm l^\prime $, we obtain
(B25)$$\begin{split}&{\sum\limits_{l^\prime = 0}^{{l_{{\max}}}} \sum\limits_{m^\prime = - l^\prime}^{{l^\prime}} \frac{{\left\langle {l^\prime |{\phi _n}(m^\prime)} \right\rangle}}{{\sqrt {{\sigma _{{l^\prime}}}}}}{Y_{l^\prime m^\prime}}(\hat{\textbf{s}};\hat{\textbf{k}})\left({- \frac{{\hat{\textbf{s}} \cdot \hat{\textbf{k}}}}{\nu} + 1} \right){\mu _t}}\\&\quad= {\mu _s}\sum\limits_{l^\prime = 0}^{{l_{{\max}}}} \sum\limits_{m^\prime = - l^\prime}^{{l^\prime}} {{\rm{g}}^{{l^\prime}}}\frac{{\left\langle {l^\prime |{\phi _n}(m^\prime)} \right\rangle}}{{\sqrt {{\sigma _{{l^\prime}}}}}}{Y_{l^\prime m^\prime}}(\hat{\textbf{s}};\hat{\textbf{k}}).\end{split}$$
By rotating the reference frame in the inverse direction, we arrive at
(B26)$$\begin{split}&{\sum\limits_{m^\prime = - {l_{{\max}}}}^{{l_{{\max}}}} \sum\limits_{l^\prime = |m^\prime |}^{{l_{{\max}}}} \frac{{\left\langle {l^\prime |{\phi _n}(m^\prime)} \right\rangle}}{{\sqrt {{\sigma _{{l^\prime}}}}}}{Y_{l^\prime m^\prime}}(\hat{\textbf{s}})\left({- \frac{{\cos \vartheta}}{\nu} + 1} \right){\mu _t}}\\&\quad= {\mu _s}\sum\limits_{m^\prime = - {l_{{\max}}}}^{{l_{{\max}}}} \sum\limits_{l^\prime = |m^\prime |}^{{l_{{\max}}}} {{\rm{g}}^{{l^\prime}}}\frac{{\left\langle {l^\prime |{\phi _n}(m^\prime)} \right\rangle}}{{\sqrt {{\sigma _{{l^\prime}}}}}}{Y_{l^\prime m^\prime}}(\hat{\textbf{s}}).\end{split}$$
By multiplying
$Y_{{lm}}^*(\hat{\textbf{s}})$ (
${-}({l_{{\max}}} - 1) \le m \le {l_{{\max}}} - 1$,
$|m| \le l \le {l_{{\max}}}$) on both sides and integrating over
$\hat{\textbf{s}} \in {{\mathbb S}^2}$, we obtain
(B27)$$\begin{split}&{\sum\limits_{l^\prime = |m|}^{{l_{{\max}}}} \left({{b_{l + 1}}(m){\delta _{l + 1,l^\prime}} + {b_l}(m){\delta _{l - 1,l^\prime}}} \right)\left\langle {l^\prime |{\phi _n}(m)} \right\rangle}\\&\quad= \frac{{{\nu _n}(m)}}{{{\mu _t}}}\left\langle {l|{\phi _n}(m)} \right\rangle .\end{split}$$
In the above equation, we wrote
$\nu = {\nu _n}(m)$. We see that
$\frac{\nu}{{{\mu _t}}} = {\nu _n}(M)/{\mu _t}$ and
$| {{\phi _n}(M)} \rangle$ are eigenvalues and eigenvectors of the following matrix-vector equation [
8,
9]:
(B28)$$B(M)\left| {{\phi _n}(M)} \right\rangle = \frac{{{\nu _n}(M)}}{{{\mu _t}}}\left| {{\phi _n}(M)} \right\rangle ,$$
where
$M = 0, \pm 1, \ldots , \pm ({l_{{\max}}} - 1)$, and matrix
$B(M) \in {{\mathbb R}^{({l_{{\max}}} - |M| + 1) \times ({l_{{\max}}} - |M| + 1)}}$ is a tridiagonal matrix whose elements are given by
(B29)$${\{B(M)\} _{{ll^\prime}}} = {b_l}(M){\delta _{l^\prime ,l - 1}} + {b_{{l^\prime}}}(M){\delta _{l^\prime ,l + 1}}$$
for
$|M| \le l,l^\prime \le {l_{{\max}}}$. We used the notation such that
$\left\langle l |B(M)| {l + 1} \right\rangle = \langle {l + 1} |B(M)| l \rangle = {b_{l + 1}}(M)$. These
${\nu _n}(M)$ are approximate eigenvalues and discretized values of the continuous spectrum of Case’s
$\nu$ [
34,
37]. We note that for each pair of
${\nu _n}(M)$ and
$\langle {l|{\phi _n}(M)} \rangle$ (
$l = |M|, \ldots ,{l_{{\max}}}$), there exists a pair of eigenvalues
${-}{\nu _n}(M)$ and eigenvectors
${(- 1)^l}\langle {l|{\phi _n}(M)} \rangle$ [
8]. For the specific intensity
$\psi (\textbf{r},\hat{\textbf{s}})$ to vanish as
$z \to \infty$, we take only
${\lfloor\left(l_{\max}-|M|+1\right)/2\rfloor}$ eigenvalues and eigenvectors such that
(B30)$${\nu _n}(M) \gt 0,\quad n = 1,2, \ldots ,\left\lfloor {\frac{{{l_{{\max}}} - |M| + 1}}{2}} \right\rfloor .$$
From the point of view of the singular eigenfunction, the method of rotated reference frames is the spherical-harmonic expansion of the singular eigenfunction [
10,
37,
38]:
(B31)$$\Phi _\nu ^m(\hat{\textbf{s}}) \approx \sum\limits_{l = |m|}^{{l_{{\max}}}} \xi _l^m(\nu){Y_{{lm}}}(\hat{\textbf{s}}).$$
Using
$\langle {{\phi _n}(M)|{\phi _n}(M)} \rangle = 1$ and
$\int_{{{\mathbb S}^2}}\mu|\Phi _\nu ^m(\hat{\textbf{s}}{)|^2} d\hat{\textbf{s}} = 2\pi {{\cal N}^m}(\nu)$ with the normalization factor
${{\cal N}^m}(\nu)$ from one-dimensional transport theory, we find
(B32)$$\xi _l^m(\nu) = \sqrt {\frac{{2\pi {\mu _t}{{\cal N}^m}(\nu)}}{{\nu {\sigma _l}}}} \left\langle {l|{\phi _n}(m)} \right\rangle .$$
Furthermore, we note that
$\left\langle {l|{\phi _{- \nu}}(M)} \right\rangle = (- {1)^l}\langle {l|{\phi _\nu}(M)} \rangle$ [
8].
The specific intensity $\psi (\textbf{r},\hat{\textbf{s}})$ is given by the superposition of eigenmodes ${\psi _\nu}(\textbf{r},\hat{\textbf{s}},\textbf{q})$ with separation constant $\nu$ as
(B33)$$\begin{split} \psi (\textbf{r},\hat{\textbf{s}}) & = \frac{1}{{{{(2\pi)}^2}}}\sum\limits_{\nu \gt 0} \int_{{{\mathbb R}^2}} {C_\nu}(\textbf{q}){\psi _\nu}(\textbf{r},\hat{\textbf{s}},\textbf{q}) {\rm{d}}\textbf{q}\\&= \frac{1}{{{{(2\pi)}^2}}}\sum\limits_{\nu \gt 0} \int_{{{\mathbb R}^2}} {C_\nu}(\textbf{q}){e^{i\textbf{q} \cdot {\boldsymbol \rho}}}\sum\limits_{l = 0}^{{l_{{\max}}}} \sum\limits_{m = - l}^l \frac{{\left\langle {l|{\phi _\nu}} \right\rangle}}{{\sqrt {{\sigma _l}}}}{{(- 1)}^m}\\&\quad\times {e^{- im{\varphi _\textbf{q}}}}d_{{mM}}^l[i\tau (\nu q)]{Y_{{lm}}}(\hat{\textbf{s}}){e^{- {\mu _t}{\hat{k}_z}(\nu q)z/\nu}} {\rm d}\textbf{q},\end{split}$$
where
${C_\nu}(\textbf{q})$ is determined later from the boundary condition. We note that
(B34)$$\begin{split} {I_{{\rm{inc}}}^{(p)}(\textbf{r},\hat{\textbf{s}})}&= {\mu _s}{e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\sum\limits_{l = 0}^{{l_{{\max}}}} \sum\limits_{m = - l}^l {\eta _{{lm}}}{e^{- im{\varphi _{{\textbf{q}_0}}}}}\\ &\quad\times \left({{R_{\mathfrak{n}}}(\hat{\textbf{s}} \cdot \hat{\textbf{z}})(- {{1)}^{l + m}} - 1} \right){Y_{{lm}}}(\hat{\textbf{s}}).\end{split}$$
Let us find
${C_\nu}(\textbf{q})$. As was done in [
11], we introduce
(B35)$$\begin{split} {{\cal B}_{{ll^\prime}}^m(\mathfrak{n})}&= \int_{{S}_ + ^2} {R_{\mathfrak{n}}}(\cos \vartheta){Y_{l^\prime m}}(\hat{\textbf{s}})Y_{{lm}}^*(\hat{\textbf{s}}) {\rm{d}}\hat{\textbf{s}}\\ &= \frac{1}{2}\sqrt {\frac{{(2l + 1)(2l^\prime + 1)(l - m)!(l^\prime - m)!}}{{(l + m)!(l^\prime + m)!}}}\\ &\quad\times \int_0^1 {R_{\mathfrak{n}}}(\mu)P_l^m(\mu)P_{{l^\prime}}^m(\mu) {\rm{d}}\mu .\end{split}$$
Note that
${\cal B}_{{ll^\prime}}^{- m}(\mathfrak{n}) = {\cal B}_{{ll^\prime}}^m(\mathfrak{n})$. Furthermore, we let
${\cal B}_{{ll^\prime}}^m(\infty)$ denote
${\cal B}_{{ll^\prime}}^m(\mathfrak{n})$, with
${R_{\mathfrak{n}}} = 1$. Let us take the Fourier transform for
${\boldsymbol \rho}$ and operate
$\int_{{\mathbb S}_ + ^2} {\rm{d}}\hat{\textbf{s}} Y_{{lm}}^*(\hat{\textbf{s}})$ on the boundary condition. By introducing
${f_{{Mn}}}(q)$ as
(B36)$${C_\nu}(\textbf{q}) = (2\pi {)^2}{f_{{Mn}}}(q)\delta (\textbf{q} - {\textbf{q}_0}),$$
we obtain
(B37)$$\begin{split} &{\sum\limits_{M = 0}^{{l_{{\max}}} - 1} \sum\limits_{n = 1}^{\lfloor({l_{{\max}}} - |M| + 1)/2\rfloor} \left[\sum\limits_{l^\prime = \max (|m|,|M|)}^{{l_{{\max}}}}\right.}\\&{\left({{\cal B}_{{ll^\prime}}^m(\infty) - {{(- 1)}^{l^\prime + m}}{\cal B}_{{ll^\prime}}^m(\mathfrak{n})} \right)\frac{{\left\langle {l^\prime |{\phi _n}(M)} \right\rangle}}{{\sqrt {{\sigma _{{l^\prime}}}}}}}\\&\quad\times {\left.\left({d_{{mM}}^{{l^\prime}}({\vartheta _{\hat{\textbf{k}}}}) + (1 - {\delta _{M0}})(- {{1)}^M}d_{m, - M}^{{l^\prime}}({\vartheta _{\hat{\textbf{k}}}})} \right)\vphantom{\sum\limits_{l^\prime = \max (|m|,|M|)}^{{l_{{\max}}}}}\right]{f_{{Mn}}}(q)}\\&\quad= {\mu _s}\sum\limits_{l^\prime = m}^{{l_{{\max}}}} {\eta _{l^\prime m}}\left({{{(- 1)}^{{l^\prime}}}{\cal B}_{{ll^\prime}}^m(\mathfrak{n}) - {{(- 1)}^m}{\cal B}_{{ll^\prime}}^m(\infty)} \right)\end{split}$$
for
$0 \le l \le {l_{{\max}}}$,
$0 \le m \le l$. Note that equations for
$m$ and
${-}m$ are the same, and hence,
${f_{- M,n}}(q) = (- {1)^M}{f_{{Mn}}}(q)$. Due to the fact that associated Legendre polynomials satisfy three-term recurrence relations, linearly independent equations are extracted from the above equations if equations with
$l = m + 1 + 2\alpha$ (
$\alpha = 0,1, \ldots , {\lfloor({l_{{\max}}} - m - 1) / 2\rfloor}$) are taken for
$m = 0,1, \ldots , {l_{{\max}}} - 1$.
Finally, we obtain
(B38)$$\psi (\textbf{r},\hat{\textbf{s}}) = {e^{i{\textbf{q}_0} \cdot {\boldsymbol \rho}}}\sum\limits_{l = 0}^{{l_{{\max}}}} \sum\limits_{m = - l}^l {(- 1)^m}{e^{- im{\varphi _{{\textbf{q}_0}}}}}{Y_{{lm}}}(\hat{\textbf{s}}){K_{{lm}}}({q_0},z),$$
where
(B39)$$\begin{split}{K_{{lm}}}({q_0},z) &= \sum\limits_{M = - ({l_{{\max}}} - 1)}^{{l_{{\max}}} - 1} \sum\limits_{n = 1}^{\lfloor({l_{{\max}}} - |M| + 1)/2\rfloor} {f_{{Mn}}}({q_0})\\&\quad\times \frac{{\left\langle {l|{\phi _n}(M)} \right\rangle}}{{\sqrt {{\sigma _l}}}}d_{{mM}}^l[i\tau ({\nu _n}(M){q_0})]{e^{- {\mu _t}{\hat{k}_z}z/{\nu _n}(M)}}.\end{split}$$
Thus,
(B40)$$\begin{split}{{J_ +}({\textbf{r}_d})}&= \int_{{\mathbb S}_ - ^2} (\cos \vartheta){I_p}({\textbf{r}_d},\hat{\textbf{s}}) {\rm{d}}\hat{\textbf{s}} + \int_{{\mathbb S}_ - ^2} (\cos \vartheta)\psi ({\textbf{r}_d},\hat{\textbf{s}}) {\rm{d}}\hat{\textbf{s}}\\ &= - {e^{i{q_0} \cdot {{\boldsymbol \rho} _d}}}{A_{{\rm{RTE}}}}({q_0}),\end{split}$$
where
(B41)$$\begin{split} {{A_{{\rm{RTE}}}}({q_0})}&= \sqrt \pi \sum\limits_{l = 0}^{{l_{{\max}}}} {{(- 1)}^l}\sqrt {2l + 1} \left({\int_0^1\mu{P_l}(\mu) {\rm{d}}\mu} \right)\\ &\quad\times \left({{\mu _s}{\eta _{l0}} + {K_l}({q_0})} \right).\end{split}$$
Here,
(B42)$$\begin{split} {{K_l}({q_0})}&= {K_{l0}}({q_0},0) = \sum\limits_{M \ge 0,n} {f_{{Mn}}}({q_0})\frac{{\left\langle {l|{\phi _n}(M)} \right\rangle}}{{\sqrt {{\sigma _l}}}}\\ &\quad\times (d_{0M}^l[i\tau ({\nu _n}(M){q_0})]\\ &\quad+ (1 - {\delta _{M0}})(- {{1)}^M}d_{0, - M}^l[i\tau ({\nu _n}(M){q_0})]).\end{split}$$
Note that
$\int_0^1\mu{P_1}(\mu) {\rm{d}}\mu = \frac{1}{3}$,
$\int_0^1\mu{P_l}(\mu) {\rm{d}}\mu = 0$ if
$l \gt 1$ is odd, and when
$l$ is even,
(B43)$$\int_0^1\mu{P_l}(\mu) {\rm{d}}\mu = \frac{{{{(- 1)}^{\frac{l}{2} + 1}}l!}}{{{2^l}(l - 1)(l + 2){{\left[{\left({\frac{l}{2}} \right)!} \right]}^2}}}.$$
Funding
Japan Society for the Promotion of Science (16K04985, 17H02081, 17H06102, 17K05572, 18H01497, 18H05240, 18K03438); HUSM Grant-in-Aid.
Acknowledgment
The Monte Carlo simulation in Section 3.A was carried out using the package MC written by Vadim A. Markel (http://whale.seas.upenn.edu/vmarkel/CODES/MC.html). The Monte Carlo eXtreme (MCX) (http://mcx.space/) was used for Monte Carlo simulations in Sections 3.B and 3.C. The solid phantom was provided by Yukari Tanikawa.
Disclosures
The authors declare no conflicts of interest.
REFERENCES
1. S. Gioux, A. Mazhar, and D. J. Cuccia, “Spatial frequency domain imaging in 2019: principles, applications, and perspectives,” J. Biomed. Opt. 24, 071613 (2019). [CrossRef]
2. S. Gioux, A. Mazhar, B. T. Lee, S. J. Lin, A. M. Tobias, D. J. Cuccia, A. Stockdale, R. Oketokoun, Y. Ashitate, E. Kelly, M. Weinmann, N. J. Durr, L. A. Moffitt, A. J. Durkin, B. J. Tromberg, and J. V. Frangioni, “First-in-human pilot study of a spatial frequency domain oxygenation imaging system,” J. Biomed. Opt. 16, 086015 (2011). [CrossRef]
3. A. Mazhar, S. A. Sharif, J. D. Cuccia, J. S. Nelson, K. M. Kelly, and A. J. Durkin, “Spatial frequency domain imaging of port wine stain biochemical composition in response to laser therapy: a pilot study,” Lasers Surg. Med. 44, 611–621 (2012). [CrossRef]
4. K. P. Nadeau, A. J. Durkin, and B. J. Tromberg, “Advanced demodulation technique for the extraction of tissue optical properties and structural orientation contrast in the spatial frequency domain,” J. Biomed. Opt. 19, 056013 (2014). [CrossRef]
5. A. Ponticorvo, D. M. Burmeister, B. Yang, B. Choi, R. J. Christy, and A. J. Durkin, “Quantitative assessment of graded burn wounds in a porcine model using spatial frequency domain imaging (SFDI) and laser speckle imaging (LSI),” Biomed. Opt. Express 5, 3467–3481 (2014). [CrossRef]
6. J. P. Angelo, S.-J. Chen, M. Ochoa, U. Sunar, S. Gioux, and X. Intes, “Review of structured light in diffuse optical imaging,” J. Biomed. Opt. 24, 071602 (2018). [CrossRef]
7. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, F. R. Ayers, and B. J. Tromberg, “Quantitation and mapping of tissue optical properties using modulated imaging,” J. Biomed. Opt. 14, 024012 (2009). [CrossRef]
8. V. A. Markel, “Modified spherical harmonics method for solving the radiative transport equation,” Waves Random Media 14, L13–L19 (2004). [CrossRef]
9. G. Panasyuk, J. C. Schotland, and V. A. Markel, “Radiative transport equation in rotated reference frames,” J. Phys. A 39, 115–137 (2006). [CrossRef]
10. M. Machida, “Singular eigenfunctions for the three-dimensional radiative transport equation,” J. Opt. Soc. Am. A 31, 67–74 (2014). [CrossRef]
11. M. Machida, G. Panasyuk, J. C. Schotland, and V. A. Markel, “The Green’s function for the radiative transport equation in the slab geometry,” J. Phys. A 43, 065402 (2010). [CrossRef]
12. A. Liemert and A. Kienle, “Spatially modulated light source obliquely incident on a semi-infinite scattering medium,” Opt. Lett. 37, 4158–4160 (2012). [CrossRef]
13. A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep. 3, 2018 (2013). [CrossRef]
14. S. Nothelfer, F. Bergmann, A. Liemert, D. Reitzle, and A. Kienle, “Spatial frequency domain imaging using an analytical model for separation of surface and volume scattering,” J. Biomed. Opt. 24, 071604 (2019). [CrossRef]
15. R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]
16. M. Machida, G. Panasyuk, J. C. Schotland, and V. A. Markel, “Diffusion approximation revisited,” J. Opt. Soc. Am. A 26, 1291–1300 (2009). [CrossRef]
17. U. Tricoli, C. M. Macdonald, A. Da Silva, and V. A. Markel, “Optimized diffusion approximation,” J. Opt. Soc. Am. A 35, 356–369 (2018). [CrossRef]
18. L. O. Svaasand, T. Spott, J. B. Fishkin, T. Pham, B. J. Tromberg, and M. W. Berns, “Reflectance measurements of layered media with diffuse photon-density waves: a potential tool for evaluating deep burns and subcutaneous lesions,” Phys. Med. Biol. 44, 801–813 (1999). [CrossRef]
19. K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, 1967).
20. K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Q. Appl. Math. 2, 164–168 (1944). [CrossRef]
21. D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963). [CrossRef]
22. M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999). [CrossRef]
23. M. Schweiger, S. R. Arridge, and I. Nissilä, “Gauss–Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol. 50, 2365–2386 (2005). [CrossRef]
24. J. J. More, B. S. Garbow, and K. E. Hillstrom, “User guide for MINPACK-1,” Argonne National Laboratory Report ANL-80-74 (1980).
25. C. K. Hayakawa, K. Karrobi, V. Pera, D. Roblyer, and V. Venugopalan, “Optical sampling depth in the spatial frequency domain,” J. Biomed. Opt. 24, 071603 (2019). [CrossRef]
26. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II. Effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal,” Appl. Opt. 42, 2915–2922 (2003). [CrossRef]
27. J. Wang, J. Lin, Y. Chen, C. G. Welle, and T. J. Pfefer, “Phantom-based evaluation of near-infrared intracranial hematoma detector performance,” J. Biomed. Opt. 24, 045001 (2019). [CrossRef]
28. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” NeuroImage 63, 921–935 (2012). [CrossRef]
29. Y. Hoshi, “Hemodynamic signals in fNIRS,” in New Horizons in Neurovascular Coupling: A Bridge Between Brain Circulation and Neural Plasticity, K. Masamoto, H. Hirase, and K. Yamada, eds., Progress in Brain Research (Elsevier, 2016), Vol. 225.
30. S. Kohno and Y. Hoshi, “Spatial distributions of hemoglobin signals from superficial layers in the forehead during a verbal-fluency task,” J. Biomed. Opt. 21, 066009 (2016). [CrossRef]
31. K. M. Case, “Elementary solutions of the transport equation and their applications,” Ann. Phys. 9, 1–23 (1960). [CrossRef]
32. N. J. McCormick and I. Kuščer, “Bi-orthogonality relations for solving half-space transport problems,” J. Math. Phys. 7, 2036–2045 (1966). [CrossRef]
33. J. R. Mika, “Neutron transport with anisotropic scattering,” Nucl. Sci. Eng. 11, 415–427 (1961). [CrossRef]
34. R. D. M. Garcia and C. E. Siewert, “On discrete spectrum calculations in radiative transfer,” J. Quant. Spectrosc. Radiat. Transfer 42, 385–394 (1989). [CrossRef]
35. R. D. M. Garcia and C. E. Siewert, “On computing the Chandrasekhar polynomials in high order and high degree,” J. Quant. Spectrosc. Radiat. Transfer 43, 201–205 (1990). [CrossRef]
36. R. D. M. Garcia and C. E. Siewert, “On the dispersion function in particle transport theory,” Z. Angew. Math. Phys. 33, 801–806 (1982). [CrossRef]
37. M. Machida, “An FN method for the radiative transport equation in three dimensions,” J. Phys. A 48, 325001 (2015). [CrossRef]
38. M. Machida, “Numerical algorithms of the radiative transport equation using rotated reference frames for optical tomography with structured illumination,” J. Quant. Spectrosc. Radiat. Transfer 234, 124–138 (2019). [CrossRef]
Open Access
Add to BibSonomy
