Inverse transport problem in fluorescence ultrasound modulated optical tomography with angularly averaged measurements

TL;DR

This paper investigates the inverse transport problem in fluorescence ultrasound modulated optical tomography with angularly averaged data, focusing on uniqueness, stability, and reconstruction algorithms, extending previous diffusion model work.

Contribution

It introduces a new analysis for the inverse transport problem in fUMOT with angularly averaged measurements, including reconstruction algorithms and numerical validation.

Findings

Proved uniqueness and stability of the inverse problem

Developed and validated reconstruction algorithms

Extended previous diffusion model analysis to transport model

Abstract

We consider an inverse transport problem in fluorescence ultrasound modulated optical tomography (fUMOT) with angularly averaged illuminations and measurements. We study the uniqueness and stability of the reconstruction of the absorption coefficient and the quantum efficiency of the fluorescent probes. Reconstruction algorithms are proposed and numerical validations are performed as well. This paper is an extension of arXiv:1804.01135, where a diffusion model for this problem was considered.

Equations

$$\begin{array}{rl}\frac{1}{c}{\partial }_{t}u(\mathbf{x},\mathbf{v},t)+\mathbf{v}\cdot \nabla u(\mathbf{x},\mathbf{v},t)+({\sigma }_{x,a}(\mathbf{x})+{\sigma }_{x,f}(\mathbf{x})+{\sigma }_{x,s}(\mathbf{x}))u(\mathbf{x},\mathbf{v},t)& \\ ={\sigma }_{x,s}(\mathbf{x}){\int }_{{\mathbb{S}}^{d-1}}p(\mathbf{v}\cdot {\mathbf{v}}^{\prime })u(\mathbf{x},{\mathbf{v}}^{\prime },t)d{\mathbf{v}}^{\prime }& \text{in }X\times {\mathbb{R}}^{+},\\ \frac{1}{c}{\partial }_{t}w(\mathbf{x},\mathbf{v},t)+\mathbf{v}\cdot \nabla w(\mathbf{x},\mathbf{v},t)+({\sigma }_{m,a}(\mathbf{x})+{\sigma }_{m,f}(\mathbf{x})+{\sigma }_{m,s}(\mathbf{x}))w(\mathbf{x},\mathbf{v},t)& \\ ={\sigma }_{m,s}(\mathbf{x}){\int }_{{\mathbb{S}}^{d-1}}p(\mathbf{v}\cdot {\mathbf{v}}^{\prime })w(\mathbf{x},{\mathbf{v}}^{\prime },t)d{\mathbf{v}}^{\prime }+S(\mathbf{x},t)& \text{in }X\times {\mathbb{R}}^{+},\\ u(\mathbf{x},\mathbf{v},t)=g(\mathbf{x},\mathbf{v},t),w(\mathbf{x},\mathbf{v},t)=0& \text{on }{\Gamma }_{-}\times {\mathbb{R}}^{+},\\ u(\mathbf{x},\mathbf{v},t)=0,w(\mathbf{x},\mathbf{v},t)=0& \text{on }X\times \{0\}.\end{array}$$

$$S(\mathbf{x},t)=\eta (\mathbf{x}){\sigma }_{x,f}(\mathbf{x}){\int }_0^t\frac{1}{\tau }{e}^{-\frac{t-s}{\tau }}\left({\int }_{{\mathbb{S}}^{d-1}}u(\mathbf{x},\mathbf{v},s)d\mathbf{v}\right)ds,$$

$$\mathbf{v}\cdot \nabla u(\mathbf{x},\mathbf{v})+({\sigma }_{x,a}+{\sigma }_{x,f}+{\sigma }_{x,s})u(\mathbf{x},\mathbf{v})$$

$$\mathbf{v}\cdot \nabla w(\mathbf{x},\mathbf{v})+({\sigma }_{m,a}+{\sigma }_{m,f}+{\sigma }_{m,s})w(\mathbf{x},\mathbf{v})$$

$$+\eta {\sigma }_{x,f}{\int }_{{\mathbb{S}}^{d-1}}u(\mathbf{x},\mathbf{v})d\mathbf{v}$$

$$u(\mathbf{x},\mathbf{v})=g(\mathbf{x},\mathbf{v})$$

$${\sigma }_{x,s}^{\epsilon }(\mathbf{x})$$

$${\sigma }_{m,s}^{\epsilon }(\mathbf{x})$$

$${\sigma }_{x,a}^{\epsilon }(\mathbf{x})$$

$${\sigma }_{m,a}^{\epsilon }(\mathbf{x})$$

$${\sigma }_{x,f}^{\epsilon }(\mathbf{x})$$

$$\mathbf{v}\cdot \nabla u_{\epsilon }(\mathbf{x},\mathbf{v})+({\sigma }_{x,a}^{\epsilon }+{\sigma }_{x,f}^{\epsilon }+{\sigma }_{x,s}^{\epsilon })u_{\epsilon }(\mathbf{x},\mathbf{v})$$

$$\mathbf{v}\cdot \nabla w_{\epsilon }(\mathbf{x},\mathbf{v})+({\sigma }_{m,a}^{\epsilon }+{\sigma }_{m,s}^{\epsilon })w_{\epsilon }(\mathbf{x},\mathbf{v})$$

$$u_{\epsilon }(\mathbf{x},\mathbf{v})=g(\mathbf{x})$$

$$\mathcal{K}f(\mathbf{x},\mathbf{v})$$

$$\mathcal{J}{u}_{\epsilon }={\int }_{{\mathbb{S}}^{d-1}}{u}_{\epsilon }(\mathbf{x},\mathbf{v})\mathbf{v}\cdot {\mathbf{n}}_{\mathbf{x}}d\mathbf{v},\mathcal{J}{w}_{\epsilon }={\int }_{{\mathbb{S}}^{d-1}}{w}_{\epsilon }(\mathbf{x},\mathbf{v})\mathbf{v}\cdot {\mathbf{n}}_{\mathbf{x}}d\mathbf{v}.$$

$$\Lambda^{\varepsilon}(\mathbf{q},\phi)=(\mathcal{J}u_{\varepsilon},\mathcal{J}w_{\varepsilon})\Big{|}_{\partial\Omega}.$$

$$A_{\epsilon }f=-({\sigma }_{x,a}^{\epsilon }+{\sigma }_{x,f}^{\epsilon }+{\sigma }_{x,s}^{\epsilon })f+{\sigma }_{x,s}^{\epsilon }\mathcal{K}f,$$

$$A_{0}f=-({\sigma }_{x,a}+{\sigma }_{x,f}+{\sigma }_{x,s})f+{\sigma }_{x,s}\mathcal{K}f;$$

$$(\mathbf{v}\cdot \nabla -A_{\epsilon })u_{\epsilon }(\mathbf{x},\mathbf{v})=0.$$

$$(-\mathbf{v}\cdot \nabla -A_0)\mathfrak{U}(\mathbf{x},\mathbf{v})=0$$

$$\mathcal{J}\mathfrak{U}(\mathbf{x}):={\int }_{{\mathbb{S}}^{d-1}}\mathbf{v}\cdot {\mathbf{n}}_{\mathbf{x}}\mathfrak{U}(\mathbf{x},\mathbf{v})d\mathbf{v}={\int }_{{\mathbb{S}}^{d-1}}\mathbf{v}\cdot {\mathbf{n}}_{\mathbf{x}}u(\mathbf{x},-\mathbf{v})d\mathbf{v}=-{\int }_{{\mathbb{S}}^{d-1}}\mathbf{v}\cdot {\mathbf{n}}_{\mathbf{x}}u(\mathbf{x},\mathbf{v})d\mathbf{v}=-\mathcal{J}u(\mathbf{x})$$

$$(A_{\epsilon }-A_0)f=\epsilon \cos (\mathbf{x}\cdot \mathbf{q}+\varphi )\left(-({\sigma }_{x,a}+{\sigma }_{x,f}+{\sigma }_{x,s})f+{\sigma }_{x,s}\mathcal{K}f\right).$$

$${\int }_{{\mathbb{S}}^{d-1}}{\int }_{\Omega }\left((A_{\epsilon }-A_0)u_{\epsilon }\right){\mathfrak{U}}_{0}d\mathbf{x}d\mathbf{v}={\int }_{{\mathbb{S}}^{d-1}}{\int }_{\Omega }\mathbf{v}\cdot \nabla (u_{\epsilon }\mathfrak{U})d\mathbf{x}d\mathbf{v}={\int }_{{\mathbb{S}}^{d-1}}{\int }_{\partial \Omega }{\mathbf{n}}_{\mathbf{x}}\cdot \mathbf{v}{u}_{\epsilon }\mathfrak{U}d\mathbf{s}d\mathbf{v}.$$

$${\int }_{{\Gamma }_{-}}{\mathbf{n}}_{\mathbf{x}}\cdot \mathbf{v}g(\mathbf{x})\mathfrak{U}d\mathbf{s}d\mathbf{v}+{\int }_{{\Gamma }_{+}}{\mathbf{n}}_{\mathbf{x}}\cdot \mathbf{v}{u}_{\epsilon }g(\mathbf{x})d\mathbf{s}d\mathbf{v}={\int }_{\partial \Omega }\left(\mathcal{J}\mathfrak{U}+\mathcal{J}{u}_{\epsilon }\right)g(\mathbf{x})d\mathbf{s}(\mathbf{x}).$$

$$u_{\epsilon }=u+\epsilon u_1+{\epsilon }^2{u}_2+\cdots .$$

$${\int }_{{\mathbb{S}}^{d-1}}{\int }_{\Omega }\cos (\mathbf{x}\cdot \mathbf{q}+\varphi )\left(-({\sigma }_{x,a}+{\sigma }_{x,f}+{\sigma }_{x,s})u+{\sigma }_{x,s}\mathcal{K}u\right)\mathfrak{U}d\mathbf{x}d\mathbf{v}.$$

$$H(\mathbf{x})$$

$$=-{\sigma }_{x,tf}{\int }_{{\mathbb{S}}^{d-1}}u(\mathbf{x},\mathbf{v})\mathfrak{U}(\mathbf{x},\mathbf{v})d\mathbf{v}+{\sigma }_{x,s}{\int }_{{\mathbb{S}}^{d-1}}\mathcal{K}u(\mathbf{x},\mathbf{v})\mathfrak{U}(\mathbf{x},\mathbf{v})d\mathbf{v}$$

$$=-{\sigma }_{x,tf}{\int }_{{\mathbb{S}}^{d-1}}u(\mathbf{x},\mathbf{v})u(\mathbf{x},-\mathbf{v})d\mathbf{v}+{\sigma }_{x,s}{\int }_{{\mathbb{S}}^{d-1}}\mathcal{K}u(\mathbf{x},\mathbf{v})u(\mathbf{x},-\mathbf{v})d\mathbf{v},$$

$$-\mathbf{v}\cdot \nabla \mathfrak{W}(\mathbf{x},\mathbf{v})+({\sigma }_{m,a}+{\sigma }_{m,s})\mathfrak{W}(\mathbf{x},\mathbf{v})$$

$$\mathfrak{W}(\mathbf{x},\mathbf{v})$$

$$-\mathbf{v}\cdot \nabla \phi (\mathbf{x},\mathbf{v})+{\sigma }_{x,tf}\phi (\mathbf{x},\mathbf{v})$$

$$\phi (\mathbf{x},\mathbf{v})$$

$$\epsilon {\int }_{\Omega }\cos (\mathbf{x}\cdot \mathbf{q}+\varphi )\left(-({\sigma }_{m,a}+{\sigma }_{m,s}){\int }_{{\mathbb{S}}^{d-1}}w(\mathbf{x},\mathbf{v})\mathfrak{W}(\mathbf{x},\mathbf{v})d\mathbf{v}+{\sigma }_{m,s}{\int }_{{\mathbb{S}}^{d-1}}\mathcal{K}w(\mathbf{x},\mathbf{v})\mathfrak{W}(\mathbf{x},\mathbf{v})d\mathbf{v}\right)d\mathbf{x}$$

$$={\int }_{{\mathbb{S}}^{d-1}}{\int }_{\Omega }\mathbf{v}\cdot \nabla (w_{\epsilon }\mathfrak{W})d\mathbf{x}d\mathbf{v}-{\int }_{\Omega }\eta {\sigma }_{x,f}^{\epsilon }(\mathcal{I}u)(\mathcal{I}\mathfrak{W})d\mathbf{x}$$

$$={\int }_{{\mathbb{S}}^{d-1}}{\int }_{\partial \Omega }{\mathbf{n}}_{\mathbf{x}}\cdot \mathbf{v}{w}_{\epsilon }\mathfrak{W}d\mathbf{s}d\mathbf{v}-{\int }_{\Omega }\eta {\sigma }_{x,f}^{\epsilon }(\mathcal{I}u)(\mathcal{I}\mathfrak{W})d\mathbf{x}.$$

$$\epsilon {\int }_{\Omega }\cos (\mathbf{x}\cdot \mathbf{q}+\varphi )\left(-{\sigma }_{x,tf}{\int }_{{\mathbb{S}}^{d-1}}u(\mathbf{x},\mathbf{v})\phi (\mathbf{x},\mathbf{v})d\mathbf{v}+{\sigma }_{x,s}{\int }_{{\mathbb{S}}^{d-1}}\mathcal{K}u(\mathbf{x},\mathbf{v})\phi (\mathbf{x},\mathbf{v})d\mathbf{v}\right)d\mathbf{x}$$

Full text

Key words. inverse transport problem, angularly averaged measurements, hybrid modality, internal data, fluorescence, ultrasound modulated optical tomography, absorption coefficient, quantum efficiency, photon currents, uniqueness, stability, reconstruction.

AMS subject classifications 2000. 35R30, 35Q60, 35J91

1 Introduction

Fluorescence optical tomography (FOT) is a popular imaging modality for biomedical and preclinical research [3, 13, 14, 16, 32]. Upon illumination by a laser pulse, fluorescent probes are exited to a metastable state and later decay to the ground state by emitting photons at a lower frequency. The emitted light and the residual excitation light are detected at the boundary for the reconstruction of the spatial concentration and lifetimes of the fluorophores.

Fluorescence ultrasound modulated optical tomography (fUMOT) is a series of FOT experiments performed under varying acoustic modulation [40, 41, 42, 30]. The acoustic modulation perturbs the optical properties of the tissue sample, allowing the measurements to provide internal information about the optical field. As the fluorescent probes have high optical contrast and tissues are acoustically homogeneous, fUMOT is expected to provide stable high contrast reconstructions with resolution comparable to the wavelength of the acoustic field. The availability of the internal data and the wellposedness of the inverse problem is generic for hybrid imaging modalities [10, 31, 26, 22, 39, 27, 24, 2].

Light propagation in tissues obeys the radiative transport equation (RTE) [23]. When the tissue environment is highly scattering, the RTE can be approximated by the diffusion equation with a suitable boundary condition [23, 4]. fUMOT in the diffusion regime has been studied in our previous work [29]. However, the diffusion approximation fails in the following two cases: when the tissue is optically thin, the characteristic length is at the same order as the transport mean free path, thus the boundary layer effect cannot be neglected; and when the scattering or the illumination source is highly anisotropic, the optical field is necessarily anisotropic near the source. In an inverse transport problem, the illuminations and measurements at the boundary can be time dependent or time independent, and angularly resolved or angularly averaged [5, 7, 12]. Time dependent measurements and angularly resolved measurements are mathematically preferable since they preserve more singularities and permit more stable and more resolved reconstruction. However, in practice, the photon transport process is too fast for accurate time dependent measurements, and angularly resolved measurements are too sensitive to noise due to possibly low particle counts in certain directions. That is, in most practical applications, time independent and angularly averaged illuminations and measurements are less expensive and more reliable [12, 5, 8].

In this paper, we study fUMOT in the radiative transport regime with time independent and angularly averaged illumination and measurements. We derive the mathematical model for fUMOT in the transport regime following the works [29, 36]. Let $u(\mathbf{x},\mathbf{v},t)$ and $w(\mathbf{x},\mathbf{v},t)$ be the excitation and emission photon densities at location $\mathbf{x}\in\Omega$, along direction $\mathbf{v}\in\mathbb{S}^{d-1}$ at time $t\in\mathbb{R}^{+}$. The governing equations of fluorescence optical tomography (FOT) are

[TABLE]

Here, $\Omega\subset\mathbb{R}^{d}(d=2,3)$ is the domain of interest, $X=\Omega\times\mathbb{S}^{d-1}$ denotes the phase space, $\Gamma_{\pm}=\{(\mathbf{x},\mathbf{v})\in\partial\Omega\times\mathbb{S}^{d-1}|\pm\mathbf{n}_{\mathbf{x}}\cdot\mathbf{v}>0\}$ are the incoming and outgoing boundary sets, $g(\mathbf{x},\mathbf{v},t)=g(\mathbf{x},\mathbf{v})\delta(t-0^{+})$ is the external excitation laser pulse, and we assume the reflection at the interface $\partial\Omega$ is negligible. $\sigma_{x,a}$ (resp. $\sigma_{m,a}$) is the intrinsic absorption coefficient of the medium at the excitation wavelength (resp. emission wavelength), $\sigma_{x,s}$ (resp. $\sigma_{m,s}$) is the intrinsic scattering coefficient of the medium at the excitation wavelength (resp. emission wavelength), and $\sigma_{x,f}$ (resp. $\sigma_{m,f}$) is the absorption coefficient of the fluorophores at the excitation wavelength (resp. emission wavelength). The emission source term $S(\mathbf{x},t)$ is proportional to the radiant energy and given by

[TABLE]

where $\eta(\mathbf{x})$ is the quantum efficiency or quantum yield of the fluorophores and $\tau$ is the fluorescence lifetime of the excited state. The integral kernel $p(\mathbf{v}\cdot\mathbf{v}^{\prime})$ is the scattering phase function, which gives the angular distribution of light intensity scattered by particle collision. With slight abuse of notation, we set $u(\mathbf{x},\mathbf{v})=\int_{0}^{\infty}u(\mathbf{x},\mathbf{v},t)dt$, $w(\mathbf{x},\mathbf{v})=\int_{0}^{\infty}w(\mathbf{x},\mathbf{v},t)dt$. Then we integrate the system (1) over time. Noticing the fact that $u(\mathbf{x},\mathbf{v},\infty)=w(\mathbf{x},\mathbf{v},\infty)=0$, we obtain a stationary RTE system for these time-integrated quantities.

[TABLE]

In practice, the coefficient $\sigma_{m,f}$ is extremely small compared to the other coefficients [37, Fig. 1.7], therefore we set it to zero hereafter. For simplicity, in what follows, we consider the isotropic illuminations only, namely $g(\mathbf{x},\mathbf{v})=g(\mathbf{x})$.

Similar to [29, 15], we consider the plane wave ultrasound modulation in the form of ${P}(\mathbf{x},t)=A\cos(\omega t)\cos(\mathbf{q}\cdot\mathbf{x}+\phi)$, where $A$ is the amplitude, $\omega$ is the frequency, $\mathbf{q}$ is the wave vector and $\phi$ is the initial phase. Under the acoustic modulation, the optical coefficients take the form [6, 11]

[TABLE]

where $\varepsilon=\frac{A\cos(\omega t)}{\rho c_{s}^{2}}\ll 1$, $\rho$ is the particle number density, and $c_{s}$ is the sound speed. Note that the time variable $t$ in $\varepsilon$ is the time on the acoustic time scale, which is approximately constant during the much faster optical process. According to [6] the quantum efficiency $\eta(\mathbf{x})$ is not modulated by the acoustic field. Combining this with the stationary RTE (3), we obtain the governing equation for fUMOT in the transport regime,

[TABLE]

where the integral operators $\mathcal{K}$ and $\mathcal{I}$ are defined as

[TABLE]

For the measurements, we record the angularly averaged boundary photon currents at both the excitation and the emission wavelengths [3, 33],

[TABLE]

For a fixed external excitation source $g$, such boundary photon currents can be measured with multiple acoustic fields with various wave vectors $\mathbf{q}$ and initial phases $\phi$. Therefore the measurement operator is

[TABLE]

The objective is to reconstruct the absorption coefficient of the fluorophores $\sigma_{x,f}(\mathbf{x})$ and the quantum efficiency $\eta(\mathbf{x})$ from the measurement operator $\Lambda^{\varepsilon}$, assuming that the unperturbed background coefficients $\sigma_{x,a}$, $\sigma_{m,a}$, $\sigma_{x,s}$ and $\sigma_{m,s}$ have been reconstructed through other imaging methods [10, 9, 34, 36].

Due to the weak coupling between $\sigma_{x,f}$ and $\eta$ in the system (5), there exists a two-step approach to simultaneously reconstructing these two coefficients. Firstly, a nonlinear inverse medium problem at the excitation wavelength is solved to recover the absorption coefficient $\sigma_{x,f}$ using internal data derived from the excitation component of the measurement operator (7). Secondly, with the knowledge of $\sigma_{x,f}$, we solve a linear inverse source problem at the emission wavelength to find the quantum efficiency $\eta$ using internal data derived from the emission component of the measurement operator (7).

For the excitation stage, we consider two scenarios: (i) For the linearized problem with some smallness assumptions, we establish existence, uniqueness, and stability results with standard transport theory; (ii) For the nonlinear problem under the assumption that $\sigma_{x,f}$ is $\alpha$-Hölder continuous and known near the boundary, we propose a proximal reconstruction method for $\sigma_{x,f}$ which algebraically depends on the internal data. The error of this reconstruction can be made arbitrarily small with a proper choice of the source, and the stability is of Lipschitz type. The key idea is to use an isotropic source that is localized around a set of points on the boundary. Under this illumination, on a line connecting two bright points on the boundary, only the ballistic part of $u(\mathbf{x},\mathbf{v})$ contributes to the leading order term of the internal data, whereas the scattering parts of $u(\mathbf{x},\mathbf{v})$ yield lower order terms. It is an analogue of the highly collimated source function in [15], where angularly resolved illuminations and measurements are allowed.

At last, we make a few comments on some relevant inverse transport problems. In the absence of acoustic modulation, inverse problems for the time independent RTE with angularly averaged measurements and illuminations are mostly open [5, 43]. The equation considered in this setting is the first equation in (3), where the sum $\sigma_{x,a}+\sigma_{x,s}+\sigma_{x,f}$ is denoted by $\sigma_{x,tf}$ [5]. When $\sigma_{x,tf}$ is unknown or $\sigma_{x,s}$ is unknown, there is no uniqueness result for the reconstruction of $\sigma_{x,tf}$ or $\sigma_{x,s}$. When only $\sigma_{x,s}$ is unknown and $\sigma_{x,tf}$ and $\sigma_{x,s}$ are small, recovering $\sigma_{x,s}$ is severely unstable [8]. In the presence of acoustic modulation, inverse problems with the time independent RTE with angularly resolved measurements are studied in [15, 6].

The rest of the paper is organized as follows. In Section 2, we extract some internal data from the measurements. A few general properties of the inverse problem are established in Section 3. In Section 4, we reconstruct $\sigma_{x,f}$ from the internal data at the excitation stage. We give results on the uniqueness and stability of $\sigma_{x,f}$ for the linearized problem, and provide an algebraic reconstruction formula for the nonlinear problem. In Section 5, assuming $\sigma_{x,f}$ has been successfully reconstructed, we recover $\eta$ from the internal data at the emission stage. The numerical experiments on synthetic data are presented in Section 6 for validation.

2 Internal data

In analogy to [6], we introduce the self-adjoint operators $A_{\varepsilon}$ and $A_{0}$ defined by

[TABLE]

then the modulated solution $u_{\varepsilon}$ satisfies

[TABLE]

We then consider the auxiliary function $\mathfrak{U}(\mathbf{x},\mathbf{v}):=u(\mathbf{x},-\mathbf{v})$, which satisfies the adjoint radiative transfer equation

[TABLE]

with $\mathfrak{U}(\mathbf{x},\mathbf{v})=g(\mathbf{x})$ on the outgoing boundary $\Gamma_{+}$. The quantity

[TABLE]

is available from the measurements. Computing $(A_{\varepsilon}-A_{0})f$ with the modulated coefficients (4), we find that

[TABLE]

Multiplying the equations (10) and (11) by $\mathfrak{U}_{0}$ and $u_{\varepsilon}$ respectively, we obtain

[TABLE]

Since the boundary illumination is isotropic, the right-hand side is equal to

[TABLE]

The right-hand side is known from the measurements by noticing that $\mathfrak{U}(\mathbf{x},\mathbf{v})=u(\mathbf{x},-\mathbf{v})$. When $\varepsilon$ is sufficiently small, we write the solution $u_{\varepsilon}$ in an asymptotic expansion

[TABLE]

Then the following quantity is known up to higher order terms in $\varepsilon$,

[TABLE]

Varying $\mathbf{q}$ and $\phi$ and performing the inverse Fourier transform, we obtain the internal data $H(\mathbf{x})$ for the excitation stage

[TABLE]

where $\sigma_{x,tf}:=\sigma_{x,t}+\sigma_{x,f}$ and $\sigma_{x,t}:=\sigma_{x,a}+\sigma_{x,s}$ denote the total absorption coefficient at the excitation wavelength with and without the fluorescence. Similarly, to compute the internal data at the emission stage, we define auxiliary functions $\mathfrak{W}$ and $\varphi$ by the equations

[TABLE]

for some strictly positive function $h(\mathbf{x})\in L^{\infty}(\partial\Omega)$, and

[TABLE]

Multipling (19) by $w_{\varepsilon}$ and (5) by $\mathfrak{W}$, we obtain

[TABLE]

Similarly, for (20) and (5), we obtain up to higher orders in $\varepsilon$,

[TABLE]

The sum of (21) and (22) gives

[TABLE]

The first term on left-hand side in (23) is known from the measurements because

[TABLE]

The second term on left-hand side is known from the boundary conditions. The third term is bounded by the Cauchy-Schwartz inequality

[TABLE]

and from Lemma 2.2 in [1],

[TABLE]

where the constant $c$ depends on $\Omega$ only. Experimentally, $\eta$ and $\sigma_{x,f}$ are usually spatially localized functions concentrated on the target cells such that $\|\eta\sigma_{x,f}\mathcal{I}\mathfrak{W}\|_{L^{2}(\Omega)}\ll 1$, hence we omit this term from (23).

Therefore, the internal data $S$ at the emission stage is

[TABLE]

Under the assumption that $\varepsilon$ is sufficiently small, the internal data $H(\mathbf{x})$ and $S(\mathbf{x})$ given by (18) and (27) for all $\mathbf{x}\in\Omega$ are available. In the diffusion regime, it is easy to check the above internal data $H$ and $S$ will be simplified to the internal data in [29]. In the following, we will recover the unknown coefficients $(\sigma_{x,f},\eta)$ from the internal data $(H,S)$ simultaneously. Since the coupling between $\sigma_{x,f}$ and $\eta$ is weak, we take a two-step reconstruction process, i.e., first reconstruct $\sigma_{x,f}$ from $H$ and then use the recovered coefficient to reconstruct the quantum efficiency $\eta$ as in [29].

3 General properties of the inverse problems

In this section, we derive some general properties of the inverse problems of reconstructing $\sigma_{x,f}$ from the internal data $H$ and reconstructing $\eta$ from $S$ in the transport equations (3). For any $1\leq p\leq\infty$, let $L^{p}(X)$ (resp. $L^{p}(\Omega)$) denote the Lebesgue space of real-valued functions whose $p$-th power are Lebesgue integrable over $X$ (resp. $\Omega$), and $\mathcal{H}_{p}^{1}(X)$ the space of $L^{p}(X)$ functions whose directional derivative along $\mathbf{v}$ belongs to $L^{p}(X)$ as well, i.e., $\mathcal{H}_{p}^{1}(X):=\{f(\mathbf{x},\mathbf{v}):f\in L^{p}(X)\text{ and }\mathbf{v}\cdot\nabla f\in L^{p}(X)\}$. We also let $L^{p}(\Gamma_{-})$ denote the space of functions that are the traces of $\mathcal{H}_{p}^{1}(X)$ functions on $\Gamma_{-}$ under the norm $\|f\|_{L^{p}(\Gamma_{-})}:=(\int_{\Gamma_{-}}|\mathbf{n}(\mathbf{x})\cdot\mathbf{v}||f|^{p}d\mathbf{v}d\mathbf{s})^{1/p}$, where $d\mathbf{s}$ is the surface measure on $\partial\Omega$. We make the flowing assumptions

1. ($\mathfrak{A}$1).

The domain $\Omega$ is convex and simply connected, and $\partial\Omega$ is $C^{2}$. 2. ($\mathfrak{A}$2).

The optical coefficients $\sigma_{x,a},\sigma_{x,s},\sigma_{m,a},\sigma_{m,s}$ are bounded by some constants $\mathsf{c}_{1}$ and $\mathsf{c}_{2}$, with

[TABLE]

The unknown coefficients $\sigma_{x,f},\eta$ belong to the admissible sets $\mathcal{A}_{\sigma}$ and $\mathcal{A}_{\eta}$, respectively, where

[TABLE]

for some constants $\mathsf{c}_{3},\mathsf{c}_{4},\mathsf{c}_{5}$ and $\mathsf{c}_{6}$. 3. ($\mathfrak{A}$3).

The source function $g(\mathbf{x})$ is strictly positive, that is, there exists a constant $\mathsf{c}_{7}$ such that $0<\mathsf{c}_{7}\leq g(\mathbf{x})$ for $\mathbf{x}\in\partial\Omega$. 4. ($\mathfrak{A}$4).

The scattering phase function $p(\mathbf{v}\cdot\mathbf{v}^{\prime})$ is strictly positive and uniformly bounded and satisfies

[TABLE]

for some constants $\mathsf{c}_{8}$ and $\mathsf{c}_{9}$.

The above assumptions permit unique solutions $u(\mathbf{x},\mathbf{v}),w(\mathbf{x},\mathbf{v})\in\mathcal{H}^{1}_{p}(X)$ to RTE (3) for any given function $g(\mathbf{x})\in L^{p}(\partial\Omega)$ from the standard transport theory in  [1]. Therefore the internal data $H$ and $S$ are well-defined for any $g(\mathbf{x})\in L^{p}(\partial\Omega)$ that satisfies the above assumptions. In the following, we show that $H$ and $S$ continuously depend on the unknown coefficients $\sigma_{x,f}$ and $\eta$ respectively.

Theorem 3.1**.**

For any $g(\mathbf{x})\in L^{p}(\partial\Omega)$, suppose the assumptions ($\mathfrak{A}$1-$\mathfrak{A}$4) hold, then the operator $\mathcal{H}:L^{\infty}(\Omega)\rightarrow L^{p/2}(\Omega)$, which maps $\sigma_{x,f}$ to the internal data $H$, is Fréchet differentiable at any $\sigma_{x,f}\in\mathcal{A}_{\sigma}$ in the direction $\delta\sigma_{x,f}\in L^{\infty}(\Omega)$ such that $\sigma_{x,f}+\delta\sigma_{x,f}\in\mathcal{A}_{\sigma}$. The derivative is given by

[TABLE]

where $v(\mathbf{x},\mathbf{v})$ satisfies

[TABLE]

Proof.

Let $\tilde{\sigma}_{x,f}=\sigma_{x,f}+\delta\sigma_{x,f}$, $\tilde{u}$ be the solution to the first equation in (3) with coefficient $\tilde{\sigma}_{x,f}$, and $\tilde{H}$ be the corresponding internal data. Then $u^{\prime}:=\tilde{u}-u$ solves the transport equation

[TABLE]

Denote the difference between $v$ and the true perturbation $u^{\prime}$ by $u^{\prime\prime}:=u^{\prime}-v$. We have that $u^{\prime\prime}$ satisfies the transport equation

[TABLE]

We now show that $\|u^{\prime\prime}\|_{L^{p}(X)}$ is of order $\|\delta\sigma_{x,f}\|^{2}_{L^{\infty}(\Omega)}$ using the standard theory of transport equations [1]. The source term $\delta\sigma_{x,f}u^{\prime}$ in (33) is in $L^{p}(X)$, therefore $u^{\prime}\in\mathcal{H}_{p}^{1}(X)$ and there exist constants $\mathfrak{c}_{1}$ and $\mathfrak{c}_{2}$ such that

[TABLE]

It follows that the source term $\delta\sigma_{x,f}u^{\prime}$ in (34) lies in $L^{p}(X)$, thus

[TABLE]

Hence $u$ is Fréchet differentiable with respect to $\sigma_{x,f}$ as a map from $L^{\infty}(\Omega)$ to $L^{p}(X)$. By the product rule, the Fréchet derivative of ${\mathcal{H}}$ with respect to $\sigma_{x,f}$ is

[TABLE]

∎

Theorem 3.2**.**

For any $g(\mathbf{x})\in L^{\infty}(\partial\Omega)$, suppose the assumptions ($\mathfrak{A}$1-$\mathfrak{A}$4) hold. Then the Fréchet derivative ${\mathcal{H}}^{\prime}[\sigma_{x,f}]:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$ is Fredholm.

Proof.

From the assumptions ($\mathfrak{A}$1-$\mathfrak{A}$4), the solution $u(\mathbf{x},\mathbf{v})$ is strictly positive, thus there exists a constant $\hat{\mathsf{c}}>0$ such that

[TABLE]

On the other hand, since $v(\mathbf{x},\mathbf{v})u(\mathbf{x},-\mathbf{v})\in\mathcal{H}_{2}^{1}(X)$ and $\mathcal{K}v(\mathbf{x},\mathbf{v})u(\mathbf{x},-\mathbf{v})\in\mathcal{H}_{2}^{1}(X)$, by the averaging lemma [17, 18], we obtain

[TABLE]

Then by the fact that the embedding from $W^{2,1/2}(\Omega)$ to $L^{2}(\Omega)$ is compact, we obtain that ${\mathcal{H}}^{\prime}[\sigma_{x,f}]:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$ is Fredholm. ∎

Theorem 3.3**.**

For any $g(\mathbf{x})\in L^{\infty}(\Omega)$, suppose the assumptions ($\mathfrak{A}$1-$\mathfrak{A}$4) hold and $\sigma_{x,f}$ is known. Then the linear operator $\mathcal{S}:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$, which maps $\eta$ to the internal data $S$, is Fredholm.

Proof.

Since $\sigma_{x,f}$ is known, $w(\mathbf{x},\mathbf{v})$ and $\varphi(\mathbf{x},\mathbf{v})$ are linear in $\eta$, hence ${\mathcal{S}}$ is a linear functional of $\eta$. Since the auxiliary function $h(\mathbf{x})$ in (19) is strictly positive, $\mathcal{I}\mathfrak{W}$ is strictly positive over $\Omega$. Thus $\sigma_{x,f}(\mathcal{I}u)(\mathcal{I}\mathfrak{W})$ is strictly positive. On the other hand, since $w(\mathbf{x},\mathbf{v})\mathfrak{W}(\mathbf{x},\mathbf{v})\in\mathcal{H}_{2}^{1}(X)$ and $u(\mathbf{x},\mathbf{v})\varphi(\mathbf{x},\mathbf{v})\in\mathcal{H}_{2}^{1}(X)$, by the averaging lemma [17, 18], we have

[TABLE]

Similarly, it is easy to verify that $\mathcal{K}w(\mathbf{x},\mathbf{v})\mathfrak{W}(\mathbf{x},\mathbf{v})\in\mathcal{H}_{2}^{1}(X)$ and $\mathcal{K}u(\mathbf{x},\mathbf{v})\psi(\mathbf{x},\mathbf{v})\in\mathcal{H}_{2}^{1}(X)$ as well, hence

[TABLE]

By the compactness of the embedding from $W^{2,1/2}(\Omega)$ into $L^{2}(\Omega)$, we obtain that ${\mathcal{S}}:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$ is Fredholm. ∎

4 Reconstruction of $\sigma_{x,f}$

In this section, we consider the reconstruction of the coefficient $\sigma_{x,f}$ from the internal data $H$ in two scenarios. We first show that the linearized inverse problem permits a unique reconstruction when the medium is optically thin and the scattering is weak. Then propose a proximal reconstruction for the nonlinear problem, which allows an arbitrary accuracy when $\sigma_{x,f}$ is Hölder continuous and is known near the boundary. Note that without these further assumptions on the medium parameters, this nonlinear inverse medium problem may not have a unique reconstruction, that is, two different $\sigma_{x,f}$’s may give the same $H$ (see Section 6).

4.1 Uniqueness and stability for linearized problem

We have the following theorem for the linearized problem.

Theorem 4.1**.**

Let $g(\mathbf{x})\in L^{\infty}(\partial\Omega)$, and suppose that the assumptions ($\mathfrak{A}$1-$\mathfrak{A}$4) hold. Let the following conditions be satisfied:

The medium is optically thin, i.e., there exists a small constant $1>\gamma>0$ such that

[TABLE] 2. 2.

The scattering is weak, i.e., there exists a small constant $1>\delta>0$ such that

[TABLE] 3. 3.

The constants $\gamma$ and $\delta$ satisfy

[TABLE]

Then the linear equation ${\mathcal{H}}^{\prime}[\sigma_{x,f}]\delta\sigma_{x,f}=0$ only permits the zero solution.

Proof.

It follows from (32) that

[TABLE]

Substituting the above equation into (31), we obtain that when ${\mathcal{H}}^{\prime}[\sigma_{x,f}](\delta\sigma_{x,f})=0$, $v(\mathbf{x},\mathbf{v})$ satisfies

[TABLE]

where the map ${\mathcal{K}_{1}}$ is defined as

[TABLE]

Define the operators $L$ and ${T}$ by

[TABLE]

Let the space $L_{\sigma}^{p}(X)$ be the space of functions with the norm

[TABLE]

By Lemma 4.1 in [38], $\|L^{-1}\|_{L^{p}_{\sigma}(X)}\leq(1-\exp(-\ell_{\Omega}\sup_{\mathbf{x}\in\Omega}\sigma_{x,tf}(\mathbf{x})))$. We also have the following estimate for $u(\mathbf{x},\mathbf{v})$ from the maximum principle and semigroup theory,

[TABLE]

Thus $\|{T}\|_{L^{p}_{\sigma}(X)}$ is bounded by

[TABLE]

From the given conditions (41), (42) and (43), we deduce that $\|L^{-1}\|_{L^{p}_{\sigma}(X)}\leq\frac{\gamma}{1+\gamma}$ and $\|{T}\|_{L^{p}_{\sigma}}\leq\delta+2\mu^{2}(1+\gamma)^{2}(1+\delta)$. Thus

[TABLE]

Therefore $Lv={T}v$ only permits $v=0$ in $L^{p}_{\sigma}(X)$, and the proof is completed by noticing $L^{p}_{\sigma}(X)$ is the same set as $L^{p}(X)$. ∎

The above local uniqueness result could be interpreted by considering the limiting case. When the scattering coefficient $\sigma_{x,s}\to 0$, the internal data $H\to-\sigma_{x,tf}\int_{\mathbb{S}^{d-1}}u(\mathbf{x},\mathbf{v})u(\mathbf{x},-\mathbf{v})d\mathbf{v}$. If the medium is optically thin or $\text{diam}(\Omega)\ll 1$, then the solution $u(\mathbf{x},\mathbf{v})$ could be well approximated by ignoring the coefficient $\sigma_{x,f}$ in (3), thus $\sigma_{x,f}$ and $\sigma_{x,tf}$ are decoupled and can be recovered directly.

The following stability estimate follows immediately from the classical stability theory of Fredholm operators [25].

Theorem 4.2**.**

Let $\mathfrak{H}$ and $\tilde{\mathfrak{H}}$ be two perturbed internal data defined in (31), and $\delta\sigma_{x,f}$ and $\delta\tilde{\sigma}_{x,f}$ be the solutions to the linearized equations

[TABLE]

where $\sigma_{x,f}$ is the background coefficient. Then under the same condition as Theorem 4.1, there exists a constant $\tilde{\mathsf{c}}=\tilde{\mathsf{c}}(\gamma,\delta)>0$ such that

[TABLE]

Remark 4.3**.**

We point out that most biological tissues are typically strongly scattering, for which $\sup_{\mathbf{x}\in\Omega}\frac{\sigma_{x,s}}{\sigma_{x,tf}}$ is close to $1$, thus the second assumption in Theorem 4.1 does not hold. For the conclusion in Theorem 4.1 to hold in this case, it requires the domain size $\ell_{\Omega}$ to be small enough. Alternatively, we introduce a proximal reconstruction method in Section 4.2, which requires neither $\sup_{\mathbf{x}\in\Omega}\frac{\sigma_{x,s}}{\sigma_{x,tf}}\ll 1$ nor $\ell_{\Omega}$ is small.

4.2 Proximal uniqueness and stability of nonlinear problem

We now show an approach to approximating the coefficient $\sigma_{x,f}$ with arbitrary accuracy when $\sigma_{x,f}$ is Hölder continuous and is known near the boundary. For preparation, we introduce the following definitions and lemma.

Definition 4.4** ($\delta$-covering).**

Let $(M,\operatorname{\mathsf{d}})$ be a metric space. The set $V$ is a $\delta$-covering of $M$ if for every $\mathbf{x}\in M$, there exists $\mathbf{y}\in V$ such that $\operatorname{\mathsf{d}}(\mathbf{y},\mathbf{x})\leq\delta$.

Definition 4.5** ($\delta$-packing).**

Let $(M,\operatorname{\mathsf{d}})$ be a metric space. The set $V$ is a $\delta$-packing of $M$ if for every $\mathbf{x}_{1}\neq\mathbf{x}_{2}\in V$, $\operatorname{\mathsf{d}}(\mathbf{x}_{1},\mathbf{x}_{2})>\delta$.

Definition 4.6**.**

Suppose $\Omega\subset\mathbb{R}^{d}$ is a convex domain and $\operatorname{\mathsf{d}}$ is a metric defined on $\mathbb{R}^{d}$. Let $V=\{{\mathbf{y}}_{i}\}_{i=1}^{n}$ be a vertex set with $\mathbf{y}_{i}\in\partial\Omega$, and $G=G(V)$ be the complete geometric graph formed from the vertices $V$. Denote the edge set of $G$ by $E$. For every $e\in E$, we define the $\theta$-tube of $e$ by

[TABLE]

We then define the $\theta$-skeleton of the graph $G$ by

[TABLE]

See Fig 1 for an illustration of the formation of $G_{\theta}(V)$.

Lemma 4.7**.**

Suppose $\Omega\subset\mathbb{R}^{d}$ is a unit ball and $\operatorname{\mathsf{d}}$ is the Euclidean metric defined on $\mathbb{R}^{d}$. Then for sufficiently small $\delta>0$, there exists a vertex set $V=\{\mathbf{y}_{i}\}_{i=1}^{n}\subset\partial\Omega$ with $n=\mathcal{O}(\delta^{1-d})$ such that the $\theta$-skeleton $G_{\theta}(V)$ generated by $V$ is a $2\delta$-covering of $\Omega$ for sufficiently small $\theta\leq\mathcal{O}(\delta^{2}/n^{2})$, i.e., for any point $\mathbf{x}\in\Omega$, there exists a point $\mathbf{y}\in G_{\theta}(V)$ such that $\operatorname{\mathsf{d}}(\mathbf{x},\mathbf{y})\leq 2\delta$.

Proof.

Given any $\delta>0$, we choose a $\delta$-packing $V$ of $\partial\Omega$ with maximal cardinality. It follows that $V$ is also a $\delta$-covering of $\partial\Omega$ and $\text{card}(V)=\mathcal{O}(\delta^{1-d})$. We claim that $G(V)$ forms a $\delta$-covering of $\Omega$. For any $\mathbf{x}\in\Omega$, we pick an arbitrary point $\mathbf{y}_{i}\in V$, and denote the other intersection of $\partial\Omega$ and the line through $\mathbf{y}_{i}$ and $\mathbf{x}$ by $\mathbf{x}^{\prime}$. Since $V$ is a $\delta$-covering of $\partial\Omega$, there exists a point $\mathbf{y}_{j}\in V$ such that $\operatorname{\mathsf{d}}(\mathbf{x}^{\prime},\mathbf{y}_{j})\leq\delta$. When $\mathbf{y}_{i}\neq\mathbf{y}_{j}$, we have that $\operatorname{\mathsf{d}}(\mathbf{x},e_{ij})\leq\operatorname{\mathsf{d}}(\mathbf{x}^{\prime},\mathbf{y}_{j})\leq\delta$, where $e_{ij}$ is the edge connecting the vertices $\mathbf{y}_{i}$ and $\mathbf{y}_{j}$ (see Fig 2). The claim is obviously true when $\mathbf{y}_{i}\neq\mathbf{y}_{j}$.

Let $\theta\leq\frac{\delta^{2}}{4n^{2}}$. Consider an edge $e_{1}\in E$. For any $e_{2}\in E$ and $e_{2}\neq e_{1}$, the length of $T_{\theta}(e_{1})\cap e_{2}$ is at most $2\theta/\sin\alpha$, where $\alpha\in(0,\frac{\pi}{2})$ is the angle between $e_{1},e_{2}$. On the other hand, since $\Omega$ is the unit ball and $|e_{1}|>\delta$ from the fact that $V$ is a $\delta$-packing of $\partial\Omega$, we must have $\sin\alpha\geq\frac{\delta}{2}$, therefore $|T_{\theta}(e_{1})\cap e_{2}|\leq\frac{4\theta}{\delta}$. Because $\text{card}(E)=\binom{n}{2}\leq n^{2}$, the total length removed from $e_{1}$ is at most $n^{2}\frac{4\theta}{\delta}\leq\delta$.

We now prove that $G_{\theta}(V)$ is a $2\delta$-covering of $\Omega$. For any $\mathbf{x}\in\Omega$, since $G(V)$ is a $\delta$-covering of $\Omega$, we can find $\mathbf{y}\in G(V)$ such that $\operatorname{\mathsf{d}}(\mathbf{x},\mathbf{y})\leq\delta$ and an edge $e\in E$ such that $\mathbf{y}\in e$. Because the total length removed from $e$ is at most $\delta$, we can find a point $\mathbf{t}\in e\cap G_{\theta}(V)$ such that $\operatorname{\mathsf{d}}(\mathbf{t},\mathbf{y})\leq\delta$ and $\operatorname{\mathsf{d}}(\mathbf{x},\mathbf{t})\leq\operatorname{\mathsf{d}}(\mathbf{x},\mathbf{y})+\operatorname{\mathsf{d}}(\mathbf{y},\mathbf{t})\leq 2\delta$. ∎

We remark that $\Omega$ is taken to be a unit ball in the above lemma only for simplicity. The proof can be easily adapted to the case when the principal curvatures of $\partial\Omega$ are bounded away from zero. In the following, we assume that the domain $\Omega$ is the unit ball and prove the global uniqueness by a constructive method. The idea is to use the fact that the quadratic term $\int_{\mathbb{S}^{d-1}}u(\mathbf{x},\mathbf{v})u(\mathbf{x},-\mathbf{v})d\mathbf{v}$ contains certain “singularities” when $g(\mathbf{x})$ is concentrated at a few points on the surface $\partial\Omega$.

Theorem 4.8**.**

Let $\Omega$ be the unit ball in $\mathbb{R}^{d}$ and $\operatorname{\mathsf{d}}$ the Euclidean metric defined on $\mathbb{R}^{d}$. Suppose the assumptions ($\mathfrak{A}$1-$\mathfrak{A}$4) hold. Let the coefficient $\sigma_{x,f}$ satisfy the following conditions:

$\sigma_{x,f}$* is $\alpha$-Hölder continuous, that is, there exists a constant $\kappa>0$, such that $\forall\mathbf{x},\mathbf{y}\in\Omega$,*

[TABLE] 2. 2.

We can decompose $\sigma_{x,f}=\sigma_{x,f}^{0}+\delta\sigma_{x,f}$, where $\sigma_{x,f}^{0}$ is the known background coefficient and the unknown $\delta\sigma_{x,f}$ has compact support in the interior subdomain $\Omega_{r}$ for some $r\in(0,1)$. Here

[TABLE]

Then for any sufficiently small $\delta>0$, we can choose an illumination source $g(\mathbf{x})\in L^{\infty}(\partial\Omega)$ such that the internal data $H$ permits a reconstruction $\widetilde{\sigma}_{x,f}$ such that $\forall\mathbf{x}\in\Omega$,

[TABLE]

Proof.

Let $\delta\ll r/5$, $n=\mathcal{O}(\delta^{1-d})$ and $\theta=\mathcal{O}(\delta^{2d})$. We construct a vertex set $V=\{\mathbf{y}_{j}\}_{j=1}^{n}\subset\partial\Omega$ whose $\theta$-skeleton $G_{\theta}(V)$ forms a $2\delta$-covering of $\Omega$ as in Lemma 4.7. Let $B(\mathbf{x},s)$ denote the ball centered at $\mathbf{x}$ with radius $s$. We consider an illumination source function $g_{h}(\mathbf{x})$ of the form

[TABLE]

where $D_{j}=B(\mathbf{y}_{j},h)\cap\partial\Omega$, $\chi_{D_{j}}$ is the characteristic function of $D_{j}$, the exponent $l=(d-1)/2$, and the parameter $h\leq\theta$ is sufficiently small such that $\{D_{j}\}_{j=1}^{n}$ are disjoint from each other. Then $G_{h}\supset G_{\theta}$ is also a $2\delta$-covering of $\Omega$ and $g_{h}\in L^{\infty}(\partial\Omega)$ for any $h>0$.

Define the operators $\mathcal{B}:L^{\infty}(\partial\Omega)\rightarrow L^{\infty}(X)$ and $\mathcal{T}:L^{\infty}(X)\rightarrow L^{\infty}(X)$ as

[TABLE]

The solution to the RTE (3) with boundary illumination source $g_{h}$ is

[TABLE]

Here $\mathcal{B}g_{h}$ is the ballistic part of the solution, $\mathcal{T}\mathcal{B}g_{h}$ is the single scattering part, and $(I-\mathcal{T})^{-1}\mathcal{T}^{2}\mathcal{B}g_{h}$ is the multiple scattering part. For each point $\mathbf{x}\in\Omega_{r-4\delta}\subset\Omega_{r/5}$, we have that

[TABLE]

It follows that $\mathcal{T}\mathcal{B}g_{h}\leq\mathcal{O}(nh^{(d-1)/2})$ and $(I-\mathcal{T})^{-1}\mathcal{T}^{2}\mathcal{B}g_{h}\leq\mathcal{O}(nh^{(d-1)/2})$. Hence for all $\mathbf{x}\in\Omega_{r}$, we can write

[TABLE]

The internal data $H_{h}$ is

[TABLE]

where $\mathfrak{g}_{h}(\mathbf{x},\mathbf{v})=g_{h}(\mathbf{x}-\tau_{-}(\mathbf{x},\mathbf{v})\mathbf{v})g_{h}(\mathbf{x}+\tau_{+}(\mathbf{x},\mathbf{v})\mathbf{v})$. Let $\mathbf{y}=\mathbf{x}+\tau_{+}(\mathbf{x},\mathbf{v})\mathbf{v}$ and $\mathbf{y}^{\prime}=\mathbf{y}-\tau_{-}(\mathbf{y},\mathbf{v})\mathbf{v}$ with $\mathbf{v}=\frac{\mathbf{y}-\mathbf{x}}{|\mathbf{y}-\mathbf{x}|}$. Utilizing the transformation

[TABLE]

with $\mathbf{n}(\mathbf{y})$ being the unit normal vector at $\mathbf{y}$, and noticing

[TABLE]

we can rewrite the formulation (60) as

[TABLE]

where $D_{ij}(\mathbf{x})=\{\mathbf{y}:\mathbf{y}\in D_{i},\mathbf{y}-\tau_{-}(\mathbf{y},\mathbf{v})\mathbf{v}\in D_{j}\;\text{with}\;\mathbf{v}=\frac{\mathbf{y}-\mathbf{x}}{|\mathbf{y}-\mathbf{x}|}\}$ and $E(\mathbf{y},\mathbf{y}^{\prime})$ is

[TABLE]

Next, since $|\mathbf{y}-\mathbf{y}_{i}|\leq h$ and $|\mathbf{y}-\tau_{-}(\mathbf{y},\mathbf{v})\mathbf{v}-\mathbf{y}_{j}|\leq h$, we take the Taylor expansion at $\mathbf{y}=\mathbf{y}_{i}$ and $\mathbf{y}^{\prime}=\mathbf{y}_{j}$ for each integral over $D_{ij}$ and obtain

[TABLE]

For an arbitrary $\mathbf{z}\in\Omega_{r-2\delta}\cap G_{h}(V)$, there exists a unique edge $e_{lk}\in E$ connecting $\mathbf{y}_{l}$ and $\mathbf{y}_{k}$ such that $\mathbf{z}=t\mathbf{y}_{k}+(1-t)\mathbf{y}_{l}$ for some $t\in(0,1)$. This means that $D_{ij}(\mathbf{z})\neq\emptyset$ if and only if $i=k,j=l$ or $i=l,j=k$. Therefore we have

[TABLE]

where $\mathfrak{B}_{h}(\mathbf{z},\mathbf{y}_{k},\mathbf{y}_{l})$ satisfies that for some constant $\mathsf{c}_{10}>0$,

[TABLE]

Here we have used $\mathbf{n}(\mathbf{y}_{k})=\mathbf{y}_{k}$ in the second equality. In general, if the principal curvatures of $\partial\Omega$ are bounded away from zero, the same lower bound in (65) still holds. Since $\mathcal{B}_{h}(\mathbf{z},\mathbf{y}_{k},\mathbf{y}_{l})\geq\mathsf{c}_{10}\delta$ for some $\mathsf{c}_{10}>0$ and $\mathcal{B}_{h}(\mathbf{z},\mathbf{y}_{k},\mathbf{y}_{l})$ is independent of $\sigma_{x,f}$, we have

[TABLE]

On the other hand, by the argument in the proof of Lemma 4.7, there exists $\mathbf{z}^{\prime}\in\Omega_{r-2\delta}^{\complement}\cap\Omega_{r-4\delta}\cap e_{lk}\cap G_{h}(V)$ (see Fig 3) such that

[TABLE]

which is known from the background coefficient $\sigma_{x,f}^{0}$. Taking the ratio between (66) and (67), we obtain

[TABLE]

Recalling that $n=\mathcal{O}(\delta^{1-d})$ and $h\leq\theta=\mathcal{O}(\delta^{2d})$, we have $\mathcal{O}\left(\frac{n^{2}h}{\delta}\right)\leq\mathcal{O}(\delta)$. Therefore for each $\mathbf{z}\in\Omega_{r-2\delta}\cap G_{h}(V)$, we can recover $\sigma_{x,tf}$ (hence $\sigma_{x,f}$) up to an error of $\mathcal{O}(\delta)$. Since $\Omega_{r-2\delta}\cap G_{h}(V)$ is a $2\delta$-covering for $\Omega_{r}$, for any $\mathbf{x}\in\Omega_{r}$, we can find a point $\mathbf{z}\in\Omega_{r-2\delta}\cap G_{h}(V)$ such that $|\mathbf{z}-\mathbf{x}|\leq 2\delta$. Using the Hölder continuity condition of $\sigma_{x,f}$, we conclude that the $L^{\infty}$ reconstruction error is bounded by $\mathcal{O}(\delta)+\kappa(2\delta)^{\alpha}=\mathcal{O}(\delta^{\alpha})$. ∎

Notice that, if the conditions in Theorem 4.8 are not satisfied, then the uniqueness of the above nonlinear case might not hold under certain circumstances. We demonstrate a numerical example which permits two distinct reconstructions for this situation in Example 6.1. In practice, the specific singular illumination source in (55) with $h\to 0$ is not possible due to resolution limitation. However, for a moderately small $h$, and a source $g_{h}$ which only concentrates at a few spots on the boundary, and when the total absorption coefficient $\sigma_{x,tf}$ is not too large, the ballistic signal still can be captured in $\int_{\mathbb{S}^{d-1}}u(\mathbf{x},\mathbf{v})u(\mathbf{x},-\mathbf{v})d\mathbf{v}$ near its $h$-skeleton. This could be used to recover the information on the $h$-skeleton approximately; see Fig 4. Although the uniqueness result of the above theorem is “proximal” and constructive, it does not rule out uniqueness for other types of illumination source.

Theorem 4.9**.**

Under the same assumptions of Theorem 4.8, let $\sigma^{1}_{x,f}$ and $\sigma^{2}_{x,f}$ be two admissible absorption coefficients of the fluorophores. Choose the illumination source $g_{h}$ constructed in Theorem 4.8 with $\delta\ll 1$, suppose $H_{1},H_{2}$ are the corresponding internal data associated with $\sigma^{1}_{x,f}$ and $\sigma^{2}_{x,f}$ respectively. If $H_{1}$ and $H_{2}$ satisfy $\left\|H_{1}/H_{2}-1\right\|_{L^{\infty}(\Omega_{r-4\delta})}<1$, then

[TABLE]

Proof.

Using the same argument as in Theorem 4.8, for any $\mathbf{z}\in\Omega_{r-2\delta}\cap G_{h}(V)$, there is a unique edge $e_{lk}\in E$ such that $\mathbf{z}\in e_{lk}$ and we can find $\mathbf{z}^{\prime}\in G_{h}(V)\cap\Omega_{r-4\delta}\cap\Omega_{r-2\delta}^{\complement}\cap e_{lk}$ such that

[TABLE]

where $\sigma^{i}_{x,tf}=\sigma_{x,a}+\sigma_{x,s}+\sigma^{i}_{x,f}$ for $i=1,2$. Taking the ratio of the above two equations and using the fact that $\sigma^{i}_{x,f}$ is known outside the subdomain $\Omega_{r}$, we obtain

[TABLE]

Since $\|H_{1}/H_{2}-1\|_{L^{\infty}(\Omega_{r-4\delta})}<1$, we obtain the error estimate for any $\mathbf{z}\in\Omega_{r-2\delta}\cap G_{h}(V)$,

[TABLE]

Again, since $(\sigma^{1}_{x,f}-\sigma^{2}_{x,f})$ is $\alpha$-Hölder continuous and $\Omega_{r-2\delta}\cap G_{h}(V)$ is a $2\delta$-covering of $\Omega_{r}$, we obtain

[TABLE]

∎

5 Reconstruction of $\eta$

Once $\sigma_{x,f}$ is reconstructed from the internal data $H$ at the excitation stage, we use the reconstructed coefficient $\sigma_{x,f}$ to reconstruct the quantum efficiency using the internal data $S$ at the emission stage. In practice, the reconstruction of $\sigma_{x,f}$ cannot be exact due to measurement noise. In the following theorem, we show that, as long as the error of the reconstructed $\sigma_{x,f}$ is controlled and a mild invertibility condition is satisfied, the error of the reconstructed $\eta$ is also controlled. This result can be understood by regarding the error in $\sigma_{x,f}$ as a perturbation of a compact operator, where the eigenvalues vary continuously with the perturbation [19].

Theorem 5.1**.**

Let $g(\mathbf{x})\in L^{\infty}(\Omega)$ and suppose that the assumptions ($\mathfrak{A}$1-$\mathfrak{A}4$) hold. Suppose $(\sigma_{x,f},\eta),(\tilde{\sigma}_{x,f},\tilde{\eta})\in\mathcal{A}_{\sigma}\times\mathcal{A}_{\eta}$ are two pairs of admissible coefficients and $\|\sigma_{x,f}-\tilde{\sigma}_{x,f}\|_{L^{\infty}(\Omega)}\leq\varepsilon^{\prime}$ is sufficiently small. Let $(u,w)$ and $(\tilde{u},\tilde{w})$ be the solutions for the coefficients $(\sigma_{x,f},\eta)$ and $(\tilde{\sigma}_{x,f},\tilde{\eta})$ respectively, and let $S$ and $\tilde{S}$ be the corresponding internal data at emission stage for $(\sigma_{x,f},\eta)$ and $(\tilde{\sigma}_{x,f},\tilde{\eta})$ respectively. Define the linear operators

[TABLE]

We then define the linear operators $\mathcal{A}_{i}:L^{2}(\Omega)\to L^{2}(\Omega)$, $1\leq i\leq 3$ as

[TABLE]

where $\mathfrak{W}$ is defined in equation (19). If zero is not an eigenvalue of $\mathcal{A}_{1}+\mathcal{A}_{2}+\mathcal{A}_{3}$, then there exists a constant $C>0$ such that

[TABLE]

Proof.

Let $\delta u(\mathbf{x},\mathbf{v})$, $\delta{w}(\mathbf{x},\mathbf{v})$ and $\delta\varphi(\mathbf{x},\mathbf{v})$ be the solutions to the following RTEs,

[TABLE]

Then we can write the solutions $\delta u,\delta\varphi,\delta w$ as

[TABLE]

where $\delta\eta=\eta-\tilde{\eta}$, $\delta\sigma_{x,f}=\sigma_{x,f}-\tilde{\sigma}_{x,f}$. We then decompose $(S-\tilde{S})$ into two parts: the first part is a Fredholm operator which acts on $\delta\eta$ and the second part is from the perturbation in $\sigma_{x,f}$,

[TABLE]

It is easy to verify that the reminder $\mathcal{R}$ has a trivial bound

[TABLE]

for some constant $C$. From the averaging lemma [17, 18], $\mathcal{A}_{1}+\mathcal{A}_{2}+\mathcal{A}_{3}$ is Fredholm. Thus if [math] is not an eigenvalue, we have the invertibility of $\mathcal{A}_{1}+\mathcal{A}_{2}+\mathcal{A}_{3}$ and

[TABLE]

for some constants $C^{\prime}$ and $C^{\prime\prime}$. ∎

6 Numerical experiments

The forward solver of RTE has been studied extensively in recent years and there are many existing numerical algorithms [35, 21, 28, 20]. In our work, we implement the forward solver by the discrete ordinate method with low order collocation scheme, where the phase space $X=\Omega\times\mathbb{S}^{d-1}$ is discretized in both spatial and angular space. In the physical space $\Omega$, we take the uniform mesh, on which the nodes are denoted by $\{\mathbf{x}_{i}\}_{i=1}^{N}$. In the angular space $\mathbb{S}^{d-1}$, we uniformly choose the angular directions $\{\mathbf{v}_{k}\}_{k=1}^{M}$ for each node $\mathbf{x}_{i}$. For a medium with weak scattering, the solution is solved quickly by the following source iteration:

[TABLE]

where $u^{T}(\mathbf{x},\mathbf{v})$ denotes the solution at the $T$-th iteration. Along each direction $\mathbf{v}_{k}$, the source iteration (80) can be computed with complete independence, hence the algorithm has a natural parallelism. Regarding the path integrals, we use the trapezoid rule for the first path integral term, which represents the ballistic contribution, and for the second path integral term, we use Simpson’s rule. The paralleled forward solver is implemented in C++ and wrapped with MATLAB’s mex interface, the source code is hosted at Github111https://github.com/lowrank/rte.

Although the Theorem 4.8 implies a constructive way to get an approximated estimate of $\sigma_{x,f}$, the singular localized sources require very fine mesh to resolve, which is not practical for numerical simulation with the discrete ordinate method. However, we still can seek for the reconstruction of $\sigma_{x,f}$ by minimizing the following objective functional:

[TABLE]

where $H^{\ast}$ is the synthetic internal data from the excitation stage and $\beta$ is the parameter of regularization. Using the linearization formula (32), we have

[TABLE]

where $\psi=\int_{\mathbb{S}^{d-1}}u(\mathbf{x},-\mathbf{v})u(\mathbf{x},\mathbf{v})d\mathbf{v}$ and $v$ is the solution to (32). The function $Q(\mathbf{x},\mathbf{v})$ is defined through

[TABLE]

We then use the quasi-Newton method (L-BFGS) to minimize the functional $J$. To simplify the evaluation process of the gradient, the adjoint state method is usually adopted. Let $q(\mathbf{x},\mathbf{v})$ be the solution to the adjoint RTE

[TABLE]

The gradient of $J$ is

[TABLE]

The quantum efficiency $\eta$ is then recovered by solving the corresponding linear inverse source problem using the reconstructed $\sigma_{x,f}$, which is

[TABLE]

where $S^{\ast}$ is the computed internal data from the emission stage and $\beta^{\prime}$ is the Tikhonov regularization parameter in case the problem is ill-posed.

In the following numerical experiments, the physical domain is $\Omega=[0,1]^{2}\subset\mathbb{R}^{2}$, and the scattering phase function is chosen as the Henyey-Greenstein’s function $p_{HG}$ in two dimension,

[TABLE]

where the constant g is the medium’s anisotropy parameter. To avoid the inverse crime, for the following numerical experiments, the synthetic data are generated on a fine discretization on both physical and angular spaces, while the inverse problems are solved on a coarse discretized phase space with roughly 600,000 unknowns. The numerical experiments are performed in MATLAB and the source code is hosted on Github 222https://github.com/lowrank/fumot-rte/..

6.1 Example 1

In this example, we demonstrate the nonuniqueness of the reconstruction of $\sigma_{x,f}$ in a medium with relatively strong scattering. Here $\sigma_{x,f}$ remains unknown on the entire domain $\Omega$. The coefficients are

[TABLE]

The illumination source is chosen as $g\equiv 1$ on the boundary and the anisotropy parameter $\texttt{g}=0.5$, the initial guess of $\sigma_{x,f}$ is set to zero. We also let the regularization parameter $\beta=0$ and assume noiseless internal data. In Fig 5, we can observe that the reconstructed image of $\sigma_{x,f}$ is completely different from the exact coefficient.

6.2 Example 2

In this example, we consider an optically thin medium where $\sigma_{x,s}$ is moderately small and $\sigma_{x,f}$ remains unknown on the entire domain $\Omega$. We set coefficients

[TABLE]

and let $\sigma_{x,f}$ be the modified Shepp-Logan phantom and $\eta$ the Derenzo phantom; see Fig 6.

The anisotropy parameter $\texttt{g}=0.5$ and the initial guess is generated randomly. The illumination source $g$ is chosen as

[TABLE]

Such source simulates the “singular” behavior in Theorem 4.8, which results with relative strong signals along the lines between the “points”. For the reconstruction, the regularization parameter is $\beta=10^{-3}$, and the internal data is polluted by a multiplicative random noise ${H}^{\ast}\leftarrow H^{\ast}(1+\tau\mathcal{U}([-1,1]))$, with $\mathcal{U}([-1,1])$ being the uniform distributed random variable and $\tau$ the noise level. The numerical reconstructions are shown in Fig 7.

After the coefficient $\sigma_{x,f}$ has been recovered, we continue to use this $\sigma_{x,f}$ to reconstruct the quantum efficiency $\eta$ from the internal data $S^{\ast}$, we also pollute the data by a multiplicative random noise $S^{\ast}\leftarrow S^{\ast}(1+\tau\mathcal{U}([-1,1]))$ with the same noise level. The Tikhonov regularization parameter is $\beta^{\prime}=10^{-8}$. The corresponding numerical reconstructions are shown in Fig 8.

7 Conclusion

In this paper, we studied the inverse problem in fluorescence ultrasound modulated optical tomography (fUMOT) in the transport regime with angularly averaged illumination and measurement. The inverse problem of interest is to recover the absorption coefficient of the fluorophores $\sigma_{x,f}$ and the quantum efficiency $\eta$.

We derived two internal functionals, $H(\mathbf{x})$ in (18) and $S(\mathbf{x})$ in (27), from the boundary measurement. Assuming knowledge of the background optical coefficients $\sigma_{x,a}$, $\sigma_{m,a}$, $\sigma_{x,s}$ and $\sigma_{m,s}$, we investigated the uniqueness and stability of the nonlinear map $\sigma_{x,f}\mapsto H$ as well as its linearization $\delta\sigma_{x,f}\mapsto\mathfrak{H}$. For the linearized map, we showed $\delta\sigma_{x,f}$ is uniquely and stably determined by $\mathfrak{H}$ for optically thin media. For the nonlinear map, we proved $\sigma_{x,f}$ can be approximately reconstructed with properly chosen illumination sources, up to an error that can be made arbitrarily small. Upon successful recovery of $\sigma_{x,f}$, we proved the quantum efficiency $\eta$ is also uniquely and stably determined by the internal functional $S$; moreover, the error in the reconstruction of $\eta$ is controllable as long as that of $\sigma_{x,f}$ is. Finally, the resulting reconstruction procedures are numerically implemented to validate the theoretical conclusions.

Acknowledgment

The research of YY was partly supported by the NSF Grant DMS-1715178, the AMS-Simons travel grant, and the start-up fund from the Michigan State University.