Fluorescence molecular imaging based on the adjoint radiative transport equation

In recent years, fluorescence imaging has received particular attention. This is due to its considerable potential for biomedical research and clinical applications [1, 2]. The quantification of fluorophore absorption and its distribution through biological tissues is of major interest. Fluorescence diffuse optical tomography (FDOT) is an imaging method that aims to reconstruct the internal distribution of fluorochromes or chromophores within biological tissues. This is based on light measurements collected at the tissue surface [1, 3–18]. Especially in small animals, this technology has facilitated monitoring of molecular activity, tumor growth, response to drug therapy, etc. In FDOT, a near-infrared excitation light source is used to measure fluorescence emission. Detection uses a CCD camera opposite the source that is rotated around the subject. However, due to the diffuse nature of light propagation in biological tissues, the image reconstruction problem is ill-posed and the images obtained are of low resolution. The reconstruction of the optical properties of the medium requires an accurate forward model for light propagation combined with an efficient inverse method.

The radiative transport equation (RTE) is considered as the gold standard for accurate prediction of light propagation through biological tissues at both the meso and the macroscale. I.e. the typical length of the scattering medium, beam diameter, etc are high compared to the wavelength of the incident beam [8–10, 13, 19]. Several studies deal with fluorescence molecular imaging based on the diffusion equation used as forward model. Here we give a non-exhaustive list of research into this topic [4, 11, 16, 18, 20–42]. The diffusion equation is an approximation relative to the RTE and has limitations in optically thin media. Equally, it fails in media with small geometries where boundary effects are dominant and where sources and detectors are not sufficiently far apart [43–46]. This presents specific problems for example in the field of small animal fluorescence imaging. Their fluorescent sources are potentially very close to detectors on the tissue surface. An image reconstruction method based on the RTE overcomes these limitations. There are very few reports in the literature concerning the inverse fluorescent source problem in fluorescence molecular imaging based on RTE as the forward model. Accordingly, this problem needs further investigation.

The image reconstruction methods for FDOT were reviewed previously [10, 47]. The methods based on the RTE can be found in [48–50] (steady-state domain) and [51–53] (frequency domain). There was also considerable work on developing other RTE-based methods for optical imaging [54–61]. The fluorescent source reconstruction is achieved by minimizing an objective function. This latter describes the discrepancy between measured light intensity and predicted numerical data on the tissue surface. The gradient of the objective function is a crucial indication of update through line minimization of the unknown fluorescent source distribution. The adjoint method is known to be efficient for inverse problems. It is used in research areas covered by shape optimization, fluid flow control, etc. The adjoint equation is derived from its primal equation using integration by parts. Gradient values can be efficiently calculated from a particular quantity of interest by solving the adjoint equation. The computation of the objective function gradient for large-scale problems using the finite differences method is extremely time consuming. A distinct advantage of using the adjoint method is that an explicit expression of the gradient is obtained that allows the gradient to be computed efficiently. This is achieved by solving an additional (adjoint) equation for the adjoint variable. The computational cost is equivalent to that of the forward model. For example, the gradient computational cost determined by the finite differences method is proportional to the number of parameters to be reconstructed. In practice, this number can be very large; as high as several thousand if the parameters vary with the spatial grid nodes. Furthermore, the adjoint method is independent from the number of parameters. Thus, it allows more accurate computing of the objective function gradient and the computational cost is much lower than for the finite differences method. Previous analyses dealing with the adjoint formulations of the RTE can be found in [62–64]. The adjoint method was applied to DOT problems based on RTE in the time domain [65–69] and the frequency domain [70–75]. In the reports cited, a forward model including partly reflected boundary conditions (Fresnel reflection) was considered in [71]. This took into account the refractive index mismatch at the air-tissue interface. The adjoint method was applied to FDOT problems based on the RTE frequency domain in [52, 53]. In [52], a unit strength isotropic source located at one point and transparent boundary conditions were considered in the forward model. Adjoint formulation of the RTE was given without a Lagrangian formulation or details. In [53], a unit strength isotropic source located at four points with partly reflected boundary conditions was considered in the forward model. The matrix form of a Lagrangian formulation was employed and the forward model was solved implicitly using a linear system. This paper proposes a reconstruction algorithm for FDOT of biological tissues. The RTE in the frequency domain is used to model light propagation. The adjoint method studied in this work provides an efficient way for solving the inverse problem of FDOT based on RTE.

The remainder of this manuscript is organized as follows. Section 2 presents the forward model in a 2D absorbing and highly-forward scattering medium, subjected to a collimated laser beam. Section 3 deals with the mathematical background of the inverse fluorescent source problem. A continuous Lagrangian formulation is used to rigorously deduce the adjoint RTE associated with partly reflected boundary conditions and an objective function gradient. The fluorescent source reconstruction algorithm is also presented. Section 4 analyzes and discusses the results to reconstruct the fluorophore absorption coefficient with simulated data in a 2D tissue-like phantom. Conclusions and future orientations are summarized in the final section.

The forward model for light propagation determines the fluorescent source distribution. This depends on an external excitation light source at the tissue surface denoted $\Upsilon^{\rm ex}$ (at the excitation wavelength of $\lambda^{\rm ex}$) and a spatial distribution of optical properties. In the frequency domain, the intensity of light source is usually modulated in the frequency range 100 MHz–1 GHz [1]. The demodulated transmitted intensities and phase shifts are measured at the tissue surface, using wavelength-dependent filters to distinguish excitation and emission signals. In this study, 2D absorbing and highly-forward scattering media like biological tissues are analyzed. The generation and propagation of fluorescent light in the medium can be accurately modeled by a set of two coupled RTEs in the frequency domain [10, 51, 52] as described below.

2.1. Introduction of some notations

To make easier the paper to read, we introduce here some notations which will be used thereafter. We assumed that the medium $\mathcal{D}$ is illuminated by a collimated laser beam of direction $\pmb{\Omega_c}$. Then, the illuminated wall of the medium is defined by:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \partial \mathcal{D}_c=\Big \{\pmb{r} \in \partial \mathcal{D},~ \pmb{\Omega_c} \cdot \pmb{n}(\pmb{r}) <0 \Big \}, \nonumber \end{align*}$$

where $\pmb{n}$ is the outward unit vector normal to the medium boundary. We also define the following sets:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Gamma^-&=\Big \{(\pmb{r},\pmb\Omega) \in \partial \mathcal{D} \times \Sigma,~ \pmb\Omega \cdot \pmb{n}(\pmb{r}) <0 \Big \},\nonumber \\ \nonumber \Gamma^+&=\Big \{(\pmb{r},\pmb\Omega) \in \partial \mathcal{D} \times \Sigma,~ \pmb\Omega \cdot \pmb{n}(\pmb{r}) >0 \Big \},\nonumber \\ \nonumber \Gamma_c^+&=\Big \{(\pmb{r},\pmb\Omega) \in \partial \mathcal{D}_c \times \Sigma,~ \pmb\Omega \cdot \pmb{n}(\pmb{r}) >0 \Big \}, \nonumber \end{align*}$$

where Σ is the interval [0, 2π].

2.2. Excitation light propagation

The external light source $\Upsilon^{\rm ex}(\pmb{r}, \omega)$, given at any location point $\pmb{r} \in \partial \mathcal{D}$ and at an angular modulation frequency $\omega \in \pmb{\mathbb{R}}$, penetrates from the outside into the medium. Part of it propagates through the medium without being deviated, while the rest is scattered in all directions. It is thus convenient to split the excitation radiance $\psi^{\rm ex}$ into two components. These are denoted $\psi_c^{\rm ex}(\pmb{r}, \omega)=\psi^{\rm ex}(\pmb{r}, \pmb{\Omega}, \omega) \delta (\pmb{\Omega}-\pmb{\Omega_c})$ for $\pmb{r} \in \mathcal{D}$ (δ is the Dirac-delta function and $\pmb{\Omega_c}$ is the direction of the collimated laser beam) and $\psi_s^{\rm ex}(\pmb{r}, \pmb{\Omega}, \omega)$ for $(\pmb{r}, \pmb\Omega) \in \mathcal{D} \times \Sigma$. They are respectively the collimated and diffuse components of excitation radiance [75, 76]. The $\psi_c^{\rm ex}(\pmb{r}, \omega)$ collimated excitation radiance is governed by the Bouguer–Beer–Lambert attenuation equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{ETR_psi_c} \Big [~\pmb\Omega_c \cdot \pmb\nabla +\Big(\frac{{\rm i}~\omega}{v^{\rm ex}}+\mu_t^{\rm ex}(\pmb{r})+\mu_a^{\rm ex \to em}(\pmb{r})\Big)\Big ]\psi_c^{\rm ex} (\pmb{r},\omega)=0~ {\rm{for}}~ \pmb{r} \in \mathcal{D}. \nonumber \end{align} \tag{ 1 }$$

The $\mu_a^{\rm ex \to em}$ coefficient is the absorption coefficient of a fluorescent source in the tissue at the excitation wavelength $\lambda^{\rm ex}$. The speed of light $v^{\rm ex}$ in the tissue, is given by the ratio $v^{\rm ex} = c/n^{\rm ex}$ of the speed of light in vacuum and the refractive index of the tissue. The boundary conditions of equation (1) are:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{BC_ETR_psi_c} \psi_c^{\rm ex}(\pmb{r},\omega)-\Upsilon^{\rm ex}(\pmb{r},\omega)=0~ {\rm{for}}~ \pmb{r} \in \partial \mathcal{D}_c~ {\rm{and}} ~ \psi_c^{\rm ex}(\pmb{r},\omega)=0 ~ {\rm{for}}~ \pmb{r} \in \partial \mathcal{D} \setminus \partial \mathcal{D}_c. \nonumber \end{align} \tag{ 2 }$$

It should be noted that (1) and (2) has an analytical solution given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{psi_c} \psi_c^{\rm ex}(\pmb{r}, \omega)=\Upsilon^{\rm ex}(\pmb{r_w},\omega)~{\rm{exp}}\Big [-\int_{\pmb{r_w}}^{\pmb{r}}\Big(\frac{{\rm i}~\omega}{v^{\rm ex}}+\mu_t^{\rm ex}(u)+\mu_a^{\rm ex \to em}(u)\Big)~{\rm d}u~\Big ]~ {\rm{for}}~ \pmb{r} \in \mathcal{D}~ {\rm{and}}~\pmb{r_w} \in \partial \mathcal{D}_c, \nonumber \end{align} \tag{ 3 }$$

where the integration is performed on the optical path in $\pmb{\Omega_c}$ direction from the point $\pmb{r_w}$ at the illuminated wall, to the point $\pmb{r}$ in the medium. The diffuse excitation radiance $\psi_s^{\rm ex}(\pmb{r}, \pmb{\Omega}, \omega)$ at location $\pmb{r} \in \mathcal{D} \subset \pmb{\mathbb{R}^2}$ in direction $\pmb{\Omega} \in \Sigma$ and angular modulation frequency $\omega \in \pmb{\mathbb{R}}$ is the solution of the RTE in the frequency domain at an excitation wavelength of $\lambda^{\rm ex}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big [~\pmb\Omega \cdot \pmb\nabla + \Big(\frac{{\rm i}~\omega}{v^{\rm ex}}+\mu_t^{\rm ex}(\pmb{r})+\mu_a^{\rm ex \to em}(\pmb{r})\Big)\Big ]\psi_s^{\rm ex} (\pmb{r},\pmb\Omega,\omega)\nonumber \\ -\mu_s^{\rm ex}(\pmb{r})~\int_{\Omega^{\prime}=2\pi}p^{\rm ex}(\pmb{\Omega^{\prime}}\cdot\pmb{\Omega})~\psi_s^{\rm ex} (\pmb{r},\pmb\Omega^{\prime},\omega)~{\rm d}\Omega^{\prime}-S_c^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)=0,\label{ETR_psi_s} \nonumber \end{align} \tag{ 4 }$$

for $(\pmb{r}, \pmb\Omega) \in \mathcal{D} \times \Sigma$ where,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{S_c} S_c^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)=~\mu_s^{\rm ex}(\pmb{r})~p^{\rm ex}(\pmb{\Omega_c}\cdot\pmb{\Omega})~\psi_c^{\rm ex}(\pmb{r},\omega), \nonumber \end{align} \tag{ 5 }$$

is an additional radiation source term to the RTE due to the diffuse part of the collimated excitation laser beam within the medium [75, 76]. The Henyey–Greenstein (H–G) phase function is the most widely-adopted scattering phase function of biomedical optics [8, 10, 76–78] and this has been used here. This function depends only on the inner product between the incident direction $\pmb{\Omega^{\prime}}$ that scattered $\pmb{\Omega}$ and the anisotropy factor $g^{\rm ex}$ (considering the excitation light). It is expressed for 2D media as:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{H_G} p^{\rm ex}(\pmb{\Omega^{\prime}} \cdot \pmb{\Omega})~ =~ \frac{1}{2\pi} \frac{1-(g^{\rm ex})^{2}}{\big(1+(g^{\rm ex})^{2}-2g^{\rm ex} ~\pmb{\Omega^{\prime}} \cdot \pmb{\Omega} \big)}. \nonumber \end{align} \tag{ 6 }$$

The boundary conditions for the diffuse excitation radiance are:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)-\rho(\Theta_{\rm sp})~\psi_s^{\rm ex}(\pmb{r},\pmb{\Omega_{\rm sp}},\omega)=0 ~ {\rm{with}}~ {\rm{cos}}~\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}= -\pmb{\Omega} \cdot \pmb{n} ~ {\rm{(specular~reflection)}},\nonumber \\ \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)-\frac{1}{\pi}\int_{\pmb{\Omega'} \cdot \pmb{n} > 0} \rho(\Theta')~ \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega'},\omega) ~\pmb{\Omega'} \cdot \pmb{n} ~{\rm d}{\Omega'}=0 ~ {\rm{with}}~ {\rm{cos}}~\Theta'= \pmb{\Omega'} \cdot \pmb{n} ~ {\rm{(diffuse~reflection)}},\label{BCETRbis} \nonumber \end{align} \tag{ 7 }$$

for $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} (\pmb{r}, \pmb\Omega) \in \Gamma^-$. The directional reflection coefficient ρ is given by Snell–Descartes laws assuming that the refractive index of the outside medium (air) is unity [10, 76, 79]. The specular reflection $\pmb{\Omega_{\rm sp}}=\pmb\Omega-2(\pmb{\Omega} \cdot \pmb{n})~\pmb{n}$ is defined as the direction from which a laser beam must hit the surface. Then, after a specular reflection it travels in the direction of $\pmb{\Omega}$ [79].

As for the excitation radiance, it is convenient to split the excitation fluence distribution into two components denoted $\Phi_c^{\rm ex}$ and $\Phi_s^{\rm ex}$. Then,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Phi_s} \Phi^{\rm ex}(\pmb{r},\omega)=\Phi_s^{\rm ex}(\pmb{r},\omega)+\Phi_c^{\rm ex}(\pmb{r},\omega)~ {\rm{for}} ~ \pmb{r} \in \mathcal{D} \nonumber \end{align} \tag{ 8 }$$

with $\Phi_s^{\rm ex}(\pmb{r}, \omega)=\int_{\Omega=2\pi}\psi_s^{\rm ex}(\pmb{r}, \pmb{\Omega}, \omega)~{\rm d}\Omega$ and $\Phi_c^{\rm ex}(\pmb{r}, \omega)=\psi_c^{\rm ex}(\pmb{r}, \omega)$.

The detector readings at the tissue surface are obtained from the exiting partial current or photon boundary flux. At each boundary point we have [10, 76]:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber Q_{\rm out}^{\rm ex}(\pmb{r},\omega)~ =~\int_{\pmb{\Omega} \cdot \pmb{n} > 0}\big [1-\rho(\Theta)\big ]~ \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega},\omega) ~\pmb{\Omega} \cdot \pmb{n} ~{\rm d}{\Omega} ~ {\rm{with}}~ {\rm{cos}}~\Theta= \pmb{\Omega} \cdot \pmb{n},\nonumber \\ \Big(+\big [1-\rho(\Theta_c)\big ]~\psi_c^{\rm ex}(\pmb{r},\omega)~\big (\pmb{\Omega_c} \cdot \pmb{n} \big) ~ {\rm{with}}~ {\rm{cos}}~\Theta_c= \pmb{\Omega_c} \cdot \pmb{n}~ {\rm{if}}~ ~ \pmb{\Omega_c} \cdot \pmb{n}(\pmb{r}) >0 \Big), \label{Q_out} \nonumber \end{align} \tag{ 9 }$$

for $\pmb{r} \in \partial \mathcal{D}$. It should be noted that the excitation reflectance $R^{\rm ex}$ is such that $R^{\rm ex}(\pmb{r}, \omega)=Q_{\rm out}^{\rm ex}(\pmb{r}, \omega)$ for $\pmb{r} \in \partial \mathcal{D}_c$.

2.3. Emission light propagation

The quantum yield by which the fluorescent source emits light in transit from an excitation state to the ground state is denoted η and τ denotes the local lifetime of the fluorescent source. The fluorescent source term is defined by [1, 34, 37]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Qmt} Q^{\rm em}(\pmb{r}, \omega)=\frac{1}{2\pi}~\frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1-{\rm i}~\omega~\tau(\pmb{r})}~\Phi^{\rm ex}(\pmb{r},\omega) ~ {\rm{for}}~ \pmb{r} \in \mathcal{D}. \nonumber \end{align} \tag{ 10 }$$

The source emits isotropic light, since all directional information is lost after excitation. Multiple fluorescent sources at a given location can be accounted for by modifying the source term. For example, the excitation fluence $\Phi^{\rm ex}$ can be presumed to be distributed among the fluorescent sources in proportion to their absorption coefficients. The equation that describes the transport of emitted light at the emission wavelength $\lambda^{\rm em}$ due to the fluorescent source in tissues is expressed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{ETRm} \Big [~\pmb\Omega \cdot \pmb\nabla +\Big(\frac{{\rm i}~\omega}{v^{\rm em}}+\mu_t^{\rm em}(\pmb{r})\Big)\Big ]\psi^{\rm em} (\pmb{r},\pmb\Omega, \omega) -\mu_s^{\rm em}(\pmb{r})~\int_{\Omega^{\prime}=2\pi}p^{\rm em}(\pmb{\Omega^{\prime}}\cdot\pmb{\Omega})~\psi^{\rm em} (\pmb{r},\pmb\Omega^{\prime}, \omega)~{\rm d}\Omega^{\prime}-Q^{\rm em}(\pmb{r},\omega)=0, \nonumber \end{align} \tag{ 11 }$$

where $\psi^{\rm em} (\pmb{r}, \pmb{\Omega}, \omega)$ denotes the emission radiance at position $\pmb{r}$ in direction $\pmb{\Omega}$ and at angular modulation frequency ω. The RTE (11) is weakly coupled to the RTE (1) via the excitation fluence distribution $\Phi^{\rm ex}$ which appears in (10). Thus, the RTE (1) has to be solved, but $\Phi^{\rm ex}$ has to be computed before the RTE (11) can be solved.

The boundary conditions for the emission radiance are:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)-\rho(\Theta_{\rm sp})~\psi^{\rm em}(\pmb{r},\pmb{\Omega_{\rm sp}},\omega)=0 ~ {\rm{with}}~ {\rm{cos}}~\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}= -\pmb{\Omega} \cdot \pmb{n} ~ {\rm{(specular~reflection)}},\nonumber \\ \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)-\frac{1}{\pi}\int_{\pmb{\Omega'} \cdot \pmb{n} > 0} \rho(\Theta')~ \psi^{\rm em}(\pmb{r},\pmb{\Omega'},\omega) ~\pmb{\Omega'} \cdot \pmb{n} ~{\rm d}{\Omega'}=0 ~ {\rm{with}}~ {\rm{cos}}~\Theta'= \pmb{\Omega'} \cdot \pmb{n} ~ {\rm{(diffuse~reflection)}}, \label{BCETRm} \nonumber \end{align} \tag{ 12 }$$

for $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} (\pmb{r}, \pmb\Omega) \in \Gamma^-$. The emission photon boundary flux is:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{Q_out_emission} Q_{\rm out}^{\rm em}(\pmb{r},\omega)~ =~\int_{\pmb{\Omega} \cdot \pmb{n} > 0 } \big [1-\rho(\Theta)\big ]~ \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega) ~\pmb{\Omega} \cdot \pmb{n} ~{\rm d}{\Omega}~ {\rm{with}}~ {\rm{cos}}~\Theta= \pmb{\Omega} \cdot \pmb{n},~ {\rm{for}}~ \pmb{r} \in \partial \mathcal{D}. \nonumber \end{align} \tag{ 13 }$$

It should be noted that the excitation reflectance $R^{\rm em}$ is such that $R^{\rm em}(\pmb{r}, \omega)=Q_{\rm out}^{\rm em}(\pmb{r}, \omega)$ for $\pmb{r} \in \partial \mathcal{D}_c$.

The inverse fluorescent source problem determines the unknown spatial source distribution $Q^{\rm em}(\pmb{r}, \omega)$. This last is proportional to the fluorophore absorption coefficient $\mu_a^{\rm ex \to em}$ in equation (10), from the measured fluorescence light intensity at wavelength $\lambda^{\rm em}$ escaping the tissue surface. The fluorophore concentration inside the tissue gives rise to its absorption coefficient as a linear function. A reconstructed map of the fluorophore absorption can be directly transposed onto a map of the biochemical probe concentration distribution. All other intrinsic tissue properties, i.e. the optical parameters $\mu_s^{\rm ex}(\pmb{r})$, $\mu_s^{\rm em}(\pmb{r})$, $\mu_a^{\rm ex}(\pmb{r})$, $\mu_a^{\rm em}(\pmb{r})$, $g^{\rm ex}$, $g^{\rm em}$, $n^{\rm ex}$, $n^{\rm em}$, are typically obtained from a previously performed reconstruction. Furthermore, the quantum yield η and the local lifetime τ of the fluorophore are both known and are usually provided by the manufacturer.

3.1. Introduction of some notations

To define more easier the state equation presented therafter, we denote by $\mathcal{R}_c^{\rm ex}(\cdot, \psi_c^{\rm ex})$, $\mathcal{R}_s^{\rm ex}(\cdot, \psi_c^{\rm ex}, \psi_s^{\rm ex})$ and $\mathcal{R}^{\rm em}(\cdot, \psi_c^{\rm ex}, \psi_s^{\rm ex}, \psi^{\rm em})$ the equations (1), (4) and (11). We also denote by $\widetilde{\mathcal{R}_c^{\rm ex}}(\cdot, \psi_c^{\rm ex})$, $\widetilde{\mathcal{R}_s^{\rm ex}}(\cdot, \psi_s^{\rm ex})$ and $\widetilde{\mathcal{R}^{\rm em}}(\cdot, \psi^{\rm em})$ the boundary conditions (2), (7) and (12). Finally, we define:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{R}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em}) = \left\{\begin{array}{@{}l@{}} \mathcal{R}_c^{\rm ex}(\theta,\psi_c^{\rm ex}) \nonumber \\ \mathcal{R}_s^{\rm ex}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex}) \nonumber \\ \mathcal{R}^{\rm em}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em}) \nonumber \\ \widetilde{\mathcal{R}_c^{\rm ex}}(\theta,\psi_c^{\rm ex}) \nonumber \\ \widetilde{\mathcal{R}_s^{\rm ex}}(\theta,\psi_s^{\rm ex}) \nonumber \\ \widetilde{\mathcal{R}^{\rm em}}(\theta,\psi^{\rm em}). \end{array}\right. \label{state_equation_first} \nonumber \end{align} \tag{ 14 }$$

We also introduce some notations related to functions spaces with their inner products. Let $L^2(\partial \mathcal{D}_c)$ be the space of the complex valued square-integrable functions on $\partial \mathcal{D}_c$. The Hermitian inner product and the norm associated with this solution space are defined by:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{boundary_inner_product} \Big \langle~f~ \Big{|} ~g~ \Big \rangle_B=\int_{\partial \mathcal{D}_c} \overline{f(\pmb{r})}~g(\pmb{r})~ {\rm d}r ~ {\rm{and}} ~ \big |\big | f \big |\big |_B^2=\Big \langle~f~ \Big{|} ~f~ \Big \rangle_B ~{\rm{for~all}}~ f,g \in L^2(\partial \mathcal{D}_c), \nonumber \end{align} \tag{ 15 }$$

where $\overline{f}$ is the complex conjugate of f.

Let $L^2(\Sigma)$ be the space of complex valued square-integrable functions on Σ associated with the Hermitian inner product:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{angular_inner_product} \Big \langle~f~ \Big{|}~ g~ \Big \rangle_A= \int_{\Omega=2\pi} f(\pmb\Omega)~g(\pmb\Omega)~{\rm d}\Omega ~ {\rm{for~all}}~ f,g \in L^2(\Sigma). \nonumber \end{align} \tag{ 16 }$$

Let $L^2(\mathcal{D})$ and $L^2(\mathcal{D} \times \Sigma)$ be the spaces of complex valued square-integrable functions on $\mathcal{D}$ and $\mathcal{D} \times \Sigma$. The two Hermitian inner products associated with these solution spaces are defined by:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{collimated_inner_products} \Big \langle~f~ \Big{|} ~g~ \Big \rangle_S=\int_{\mathcal{D}} \overline{f(\pmb{r})}~g(\pmb{r})~ {\rm d}r~ ; ~ \Big \langle~k~ \Big{|}~ h~ \Big \rangle_{SA}= \int_{\mathcal{D}} \int_{\Omega=2\pi} \overline{k(\pmb{r},\pmb\Omega)}~h(\pmb{r},\pmb\Omega)~{\rm d}\Omega~{\rm d}r, \nonumber \end{align} \tag{ 17 }$$

for all $f, g \in L^2(\mathcal{D})$ and $k, h \in L^2(\mathcal{D} \times \Sigma)$.

Let $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} L^2(\Gamma^-)$ be the space of complex valued square-integrable functions on $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} \Gamma^-$ associated with the Hermitian inner product:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \newcommand{\B}{{\bf B}} \displaystyle \label{boundary_inner_product_B_} \Big \langle~f~ \Big{|} ~g~ \Big \rangle_{BA^-}=\int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{f(\pmb{r},\pmb\Omega)}~g(\pmb{r},\pmb\Omega)~{\rm d}\Omega~ {\rm d}r~ {\rm{for~all}}~ f,g \in L^2(\Gamma^-). \nonumber \end{align} \tag{ 18 }$$

3.2. The objective function and observation equation

The spatial distribution of the fluorophore absorption coefficient $\mu_a^{\rm ex \to em}$ is reconstructed by applying a nonlinear optimization technique to an objective function J that is an explicit function of $\theta=\mu_a^{\rm ex \to em}$. The real-value objective function describes the discrepancy between the measured light intensity, $d_{\rm obs}^{\rm em}(\pmb{r})$ and the predicted numerical data, $R^{\rm em}(\pmb{r})$ (given from equation (13)) for positions $\pmb{r} \in \partial \mathcal{D}_c$:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{function_objective} J(\theta)=\frac{1}{2} \big |\big | R^{\rm em}(\theta) -d_{\rm obs}^{\rm em} \big |\big |_B^2 \nonumber \end{align} \tag{ 19 }$$

where $R^{\rm em}(\theta)$ is the observation equation defined as:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{observation_equation} (H \psi_{\theta}^{\rm em})(\pmb{r},\omega)= R^{\rm em}(\theta)(\pmb{r},\omega)=\int_{\pmb{\Omega} \cdot \pmb{n} > 0 } \big [1-\rho(\Theta)\big ]~ \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega) ~\pmb{\Omega} \cdot \pmb{n} ~{\rm d}{\Omega}~ {\rm{with}}~ {\rm{cos}}~\Theta= \pmb{\Omega} \cdot \pmb{n},~ {\rm{for}}~ \pmb{r} \in \partial \mathcal{D}_c. \nonumber \end{align} \tag{ 20 }$$

The reconstruction algorithm consists of minimizing J when the state equation is satisfied [72]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{state_equation} \mathcal{R}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em}) = 0. \nonumber \end{align} \tag{ 21 }$$

The nonlinear optimization algorithm chosen in this work requires knowledge of the objective function gradient with respect to unknowns. To compute this gradient, the adjoint method is introduced as described in the following section.

3.3. The Lagrangian and adjoint method

The Lagrangian is written in the L² space as [80, 81]:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber &\qquad \qquad \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em}, \widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})=\frac{1}{2} \big |\big | (H \psi^{\rm em}) -d_{\rm obs}^{\rm em} \big |\big |_B^2\nonumber \\ &+\Big \langle\phi_c^{\rm ex} \Big{|} \mathcal{R}_c^{\rm ex} \Big \rangle_S + \Big \langle \phi_s^{\rm ex} \Big{|} \mathcal{R}_s^{\rm ex} \Big \rangle_{SA}+\Big \langle\phi^{\rm em} \Big{|} \mathcal{R}^{\rm em} \Big \rangle_{SA}+\Big \langle\widetilde{\phi_c^{\rm ex}} \Big{|} \widetilde{\mathcal{R}_c^{\rm ex}} \Big \rangle_B+ \Big \langle\widetilde{\phi_s^{\rm ex}} \Big{|} \widetilde{\mathcal{R}_s^{\rm ex}} \Big\rangle_{BA^-}+\Big \langle\widetilde{\phi^{\rm em}} \Big{|} \widetilde{\mathcal{R}^{\rm em}} \Big \rangle_{BA^-},\label{Lagrangian} \nonumber \end{align} \tag{ 22 }$$

where the Lagrangian multipliers are: $\phi_c^{\rm ex}=\phi_c^{\rm ex}(\pmb{r}, \omega)$ (with $\pmb{r} \in \mathcal{D}$), $\phi_s^{\rm ex}=\phi_s^{\rm ex}(\pmb{r}, \pmb\Omega, \omega)$ and $\phi^{\rm em}=\phi^{\rm em}(\pmb{r}, \pmb\Omega, \omega)$ (with $(\pmb{r}, \pmb\Omega) \in \mathcal{D} \times \Sigma$), $\widetilde{\phi_c^{\rm ex}}=\widetilde{\phi_c^{\rm ex}}(\pmb{r}, \omega)$ (with $\pmb{r} \in \partial\mathcal{D}_c$), $\widetilde{\phi_s^{\rm ex}}=\widetilde{\phi_s^{\rm ex}}(\pmb{r}, \pmb\Omega, \omega)$ and $\widetilde{\phi^{\rm em}}=\widetilde{\phi^{\rm em}}(\pmb{r}, \pmb\Omega, \omega)$ (with $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} (\pmb{r}, \pmb\Omega) \in \Gamma^-$). They are complex functions that represent the adjoint variables associated to $\psi_c^{\rm ex}$, $\psi_s^{\rm ex}$, $\psi^{\rm em}$ in the medium and at the medium boundary. The fundamental remark (which is trivial in appearance) is that if ($\psi_c^{\rm ex}, \psi_s^{\rm ex}, \psi^{\rm em}$) is the solution of the state equation (21) for the θ parameter, then we have the identity:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{trivial_identity} \mathcal{L}(\theta,\psi_c^{\rm ex}(\theta),\psi_s^{\rm ex}(\theta),\psi^{\rm em}(\theta),\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})=J(\theta), ~{\rm{for~all}}~\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}}. \nonumber \end{align} \tag{ 23 }$$

By deriving this equation it yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber &J'(\theta)~\delta \theta=\frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \theta} \delta \theta\nonumber \\ \nonumber &~~~~~~~~~~~~~~~+\frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \psi_c^{\rm ex}} \frac{\partial \psi_c^{\rm ex}(\theta)}{\partial \theta} \delta \theta\nonumber \\ \nonumber &~~~~~~~~~~~~~~~+\frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \psi_s^{\rm ex}} \frac{\partial \psi_s^{\rm ex}(\theta)}{\partial \theta} \delta \theta\nonumber \\ &~~~~~~~~~~~~~~~~~+\frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \psi^{\rm em}} \frac{\partial \psi^{\rm em}(\theta)}{\partial \theta} \delta \theta. \label{deriving_trivial_identity} \nonumber \end{align} \tag{ 24 }$$

We denote the following independant quantities by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{independant_quantities} \delta \psi_c^{\rm ex}=\frac{\partial \psi_c^{\rm ex}(\theta)}{\partial \theta} \delta \theta,~ \delta \psi_s^{\rm ex}=\frac{\partial \psi_s^{\rm ex}(\theta)}{\partial \theta} \delta \theta ~ {\rm{and}} ~ \delta \psi^{\rm em}=\frac{\partial \psi^{\rm em}(\theta)}{\partial \theta} \delta \theta. \nonumber \end{align} \tag{ 25 }$$

Then, the adjoint variables are solutions to the following equation [80, 81]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber &\frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex}\nonumber \\ \nonumber &+\frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex}\nonumber \\ &+\frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \psi^{\rm em}} \delta \psi^{\rm em}=0,\label{abstract_adjoint_state} \nonumber \end{align} \tag{ 26 }$$

and equation (24) is reduced to:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{deriving_trivial_identity_bis} J'(\theta)~\delta \theta=\Big \langle\pmb\nabla J (\theta) \Big {|} \delta \theta \Big \rangle_S= \frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em},\widetilde{\phi_c^{\rm ex}},\widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}})}{\partial \theta} \delta \theta. \nonumber \end{align} \tag{ 27 }$$

Using (22) and (26) we deduce that:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle\phi_c^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_c^{\rm ex}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_S +\Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{SA} +\Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA} +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{SA} \nonumber \\ \nonumber +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{SA} +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA} +\Big \langle(H \psi^{\rm em})-d_{\rm obs}^{\rm em} \Big{|} (H \delta\psi^{\rm em}) \Big \rangle_B\nonumber \\ +\Big \langle\widetilde{\phi_c^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_c^{\rm ex}} }{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_B +\Big \langle\widetilde{\phi_s^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_s^{\rm ex}} }{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big\rangle_{BA^-} +\Big \langle\widetilde{\phi^{\rm em}} \Big{|} \frac {\partial \widetilde{\mathcal{R}^{\rm em}} }{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{BA^-}=0. \label{abstract_adjoint_state_2} \nonumber \end{align} \tag{ 28 }$$

Using (20), we change $(H \psi^{\rm em})$ by $R^{\rm em}$ in the right-hand side term in the second line of (28). As equation (28) has to be satisfied for all sensitivity directions $\delta \psi_c^{\rm ex}$, $\delta \psi_s^{\rm ex}$ and $\delta \psi^{\rm em}$, then it leads to the following set of equations (for each sensitivity directions):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{SA} +\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \delta\psi^{\rm em}) \Big \rangle_B +\Big \langle\widetilde{\phi^{\rm em}} \Big{|} \frac {\partial \widetilde{\mathcal{R}^{\rm em}} }{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{BA^-}=0,\nonumber \\ \nonumber \Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{SA} +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{SA} \nonumber +\Big \langle\widetilde{\phi_s^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_s^{\rm ex}} }{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big\rangle_{BA^-}=0,\nonumber \\ \Big \langle\phi_c^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_c^{\rm ex}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_S +\Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA} +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA} +\Big \langle\widetilde{\phi_c^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_c^{\rm ex}} }{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_B=0. \label{differentiation_1_bis} \nonumber \end{align} \tag{ 29 }$$

The last terms of each equation of (29) are related to the boundary conditions. Then, (29) leads also to:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{SA} +\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \delta\psi^{\rm em}) \Big \rangle_B=0,\nonumber \\ \nonumber \Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{SA} \nonumber +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{SA}=0,\nonumber \\ \Big \langle\phi_c^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_c^{\rm ex}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_S +\Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA} \nonumber +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA}=0, \nonumber \\ \Big \langle\widetilde{\phi_c^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_c^{\rm ex}} }{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_B=0;~ \Big \langle\widetilde{\phi_s^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_s^{\rm ex}} }{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big\rangle_{BA^-}=0;~ \Big \langle\widetilde{\phi^{\rm em}} \Big{|} \frac {\partial \widetilde{\mathcal{R}^{\rm em}} }{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{BA^-}=0. \label{differentiation_1} \nonumber \end{align} \tag{ 30 }$$

We denote $\mathcal{A}^*$ the adjoint operator of $\mathcal{A}$. Using the definition of the adjoint operator, the equations of (30) lead to:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi^{\rm em}}\Big)^* \phi^{\rm em} \Big{|} \delta \psi^{\rm em} \Big \rangle_{SA} +\Big \langle H^* \Big(R^{\rm em}-d_{\rm obs}^{\rm em} \Big) \Big{|} \delta \psi^{\rm em} \Big \rangle_B=0,\nonumber \\ \nonumber \Big \langle \Big(\frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_s^{\rm ex}} \Big)^* \phi_s^{\rm ex} \Big{|} \delta \psi_s^{\rm ex} \Big \rangle_{SA} +\Big \langle \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}} \Big)^* \phi^{\rm em} \Big{|} \delta \psi_s^{\rm ex} \Big \rangle_{SA}=0,\nonumber \\ \nonumber \Big \langle \Big(\frac {\partial \mathcal{R}_c^{\rm ex}}{\partial \psi_c^{\rm ex}} \Big)^* \phi_c^{\rm ex} \Big{|} \delta \psi_c^{\rm ex} \Big \rangle_S +\Big \langle \Big(\frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}} \Big)^{\overline{*}} \phi_s^{\rm ex} \Big{|} \delta \psi_c^{\rm ex} \Big \rangle_S +\Big \langle \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}}\Big)^{\overline{*}} \phi^{\rm em} \Big{|} \delta \psi_c^{\rm ex} \Big \rangle_S=0, \nonumber \\ \Big \langle \Big(\frac {\partial \widetilde{\mathcal{R}_c^{\rm ex}} }{\partial \psi_c^{\rm ex}} \Big)^* \widetilde{\phi_c^{\rm ex}} \Big{|} \delta \psi_c^{\rm ex} \Big \rangle_B=0;~ \Big \langle \Big (\frac {\partial \widetilde{\mathcal{R}_s^{\rm ex}} }{\partial \psi_s^{\rm ex}} \Big)^* \widetilde{\phi_s^{\rm ex}} \Big{|} \delta \psi_s^{\rm ex} \Big\rangle_{BA^-}=0;~ \Big \langle \Big(\frac {\partial \widetilde{\mathcal{R}^{\rm em}} }{\partial \psi^{\rm em}}\Big)^* \widetilde{\phi^{\rm em}} \Big{|} \delta \psi^{\rm em} \Big \rangle_{BA^-}=0, \label{differentiation_2} \nonumber \end{align} \tag{ 31 }$$

where $\mathcal{A}^{\overline{*}}=2 \pi ~\mathcal{A}^*$. As the equations of (31) have to be satisfied for all sensitivity directions $\delta \psi_c^{\rm ex}$, $\delta \psi_s^{\rm ex}$ and $\delta \psi^{\rm em}$, the adjoint variables must be solutions to the following set of equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi^{\rm em}}\Big)^* \phi^{\rm em} \Big{|} \delta \psi^{\rm em} \Big \rangle_{SA} +\Big \langle H^* \Big(R^{\rm em}-d_{\rm obs}^{\rm em} \Big) \Big{|} \delta \psi^{\rm em} \Big \rangle_B=0,\nonumber \\ \nonumber \Big(\frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_s^{\rm ex}} \Big)^* \phi_s^{\rm ex}+\Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}} \Big)^* \phi^{\rm em} =0,\nonumber \\ \nonumber \Big(\frac {\partial \mathcal{R}_c^{\rm ex}}{\partial \psi_c^{\rm ex}} \Big)^* \phi_c^{\rm ex} + \Big(\frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}} \Big)^{\overline{*}} \phi_s^{\rm ex} +\Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}}\Big)^{\overline{*}} \phi^{\rm em}=0, \nonumber \\ \Big(\frac {\partial \widetilde{\mathcal{R}_c^{\rm ex}} }{\partial \psi_c^{\rm ex}} \Big)^* \widetilde{\phi_c^{\rm ex}} =0;~ \Big(\frac {\partial \widetilde{\mathcal{R}_s^{\rm ex}} }{\partial \psi_s^{\rm ex}} \Big)^* \widetilde{\phi_s^{\rm ex}} =0;~ \Big(\frac {\partial \widetilde{\mathcal{R}^{\rm em}} }{\partial \psi^{\rm em}}\Big)^* \widetilde{\phi^{\rm em}} =0. \label{differentiation_3} \nonumber \end{align} \tag{ 32 }$$

Using appendices A and B and replacing $\mathcal{R}_c^{\rm ex}$, $\mathcal{R}_s^{\rm ex}$ and $\mathcal{R}^{\rm em}$ defined by (1), (4) and (11) in equation (32), we obtain the following adjoint equations model:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big [~\pmb\Omega \cdot \pmb\nabla +\Big(\frac{-{\rm i}~\omega}{v^{\rm em}}+u_t^{\rm em}(\pmb{r})\Big)\Big ]\phi^{\rm em} (\pmb{r},-\pmb\Omega, \omega) =\mu_s^{\rm em}(\pmb{r})~\int_{\Omega^{\prime}=2\pi}p^{\rm em}(\pmb{\Omega^{\prime}}\cdot (-\pmb{\Omega}))~\phi^{\rm em} (\pmb{r},\pmb\Omega^{\prime}, \omega)~{\rm d}\Omega^{\prime},\nonumber \\ \nonumber \Big [~\pmb\Omega \cdot \pmb\nabla + \Big(\frac{-{\rm i}~\omega}{v^{\rm ex}}+\mu_t^{\rm ex}(\pmb{r})+\mu_a^{\rm ex \to em}(\pmb{r})\Big)\Big ]\phi_s^{\rm ex} (\pmb{r},-\pmb\Omega,\omega)\nonumber \\ \nonumber =\mu_s^{\rm ex}(\pmb{r})~\int_{\Omega^{\prime}=2\pi}p^{\rm ex}(\pmb{\Omega^{\prime}}\cdot (-\pmb{\Omega}))~\phi_s^{\rm ex} (\pmb{r},\pmb\Omega^{\prime},\omega)~{\rm d}\Omega^{\prime}+ \frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1+{\rm i}~\omega~\tau(\pmb{r})}~\phi^{\rm em}(\pmb{r},-\pmb{\Omega},\omega),\nonumber \\ \nonumber \Big [~\pmb\Omega_c \cdot \pmb\nabla +\Big(\frac{-{\rm i}~\omega}{v^{\rm ex}}+\mu_t^{\rm ex}(\pmb{r})+\mu_a^{\rm ex \to em}(\pmb{r})\Big)\Big ]\phi_c^{\rm ex} (\pmb{r},\omega) \nonumber \\ =\mu_s^{\rm ex}(\pmb{r})~\int_{\Omega^{\prime}=2\pi}p^{\rm ex}(\pmb{\Omega^{\prime}}\cdot (\pmb{-\Omega_c}))~\phi_s^{\rm ex} (\pmb{r},\pmb\Omega^{\prime},\omega)~{\rm d}\Omega^{\prime}+\frac{1}{2\pi}~\frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1+{\rm i}~\omega~\tau(\pmb{r})}~\int_{\Omega^{\prime}=2\pi}\phi^{\rm em}(\pmb{r},\pmb{\Omega^{\prime}},\omega)~{\rm d}\Omega^{\prime}. \label{adjoint_equations_model_moins} \nonumber \end{align} \tag{ 33 }$$

The directions $\pmb{\Omega}$ and $\pmb{\Omega_c}$ were changed to $-\pmb{\Omega}$ and $-\pmb{\Omega_c}$ for convenience. The boundary conditions of the set of equation (33) were derived as shown in appendices B and C and are summarized below. If the medium boundary is assumed to be transparent, the boundary conditions of the first two equations of system (33) are given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{BC_adjoint_equations_phi_m_transparent_boundary} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)= 0 ~{\rm{for}}~ (\pmb{r},\pmb\Omega) \in \Gamma^+ \setminus \Gamma_c^+~ {\rm{and}}~ ~ \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)=- (R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma_c^+, \nonumber \end{align} \tag{ 34 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{BC_adjoint_equations_phi_c_x_phi_s_x_transparent_boundary} \phi_s^{\rm ex}(\pmb{r},\pmb\Omega,\omega)=0 ~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma^+. \nonumber \end{align} \tag{ 35 }$$

Otherwise, the medium boundary is semi-transparent with boundary conditions given by (7) and (12). Then, the adjoint boundary conditions in the case of diffuse reflection of the first two equations of system (33) are:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \newcommand{\B}{{\bf B}} \displaystyle \label{BC_adjoint_equations_phi_m} \phi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)=\frac{1}{\pi}\int_{\pmb{\Omega'} \cdot \pmb{n} < 0} \rho(\Theta')~ \phi^{\rm em}(\pmb{r},\pmb{\Omega'},\omega) ~\Big{|} \pmb{\Omega'} \cdot \pmb{n} \Big{|} ~{\rm d}{\Omega'} ~ {\rm{for}}~ (\pmb{r},\pmb\Omega) \in \Gamma^+ \setminus \Gamma_c^+, \nonumber \end{align} \tag{ 36 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \newcommand{\B}{{\bf B}} \displaystyle \label{BC_adjoint_equations_phi_m_bis} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)=& -\big [1-\rho(\Theta)\big ]~(R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega) \nonumber \\ & +\frac{1}{\pi}\int_{\pmb\Omega' \cdot \pmb n <0} \rho(\Theta')~\phi^{\rm em}(\pmb{r},\pmb\Omega',\omega)~\Big{|}\pmb\Omega' \cdot \pmb n \Big{|}~{\rm d}\Omega'~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma_c^+, \nonumber \end{align} \tag{ 37 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \newcommand{\B}{{\bf B}} \displaystyle \label{BC_adjoint_equations_phi_s} \phi_s^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)=\frac{1}{\pi}\int_{\pmb{\Omega'} \cdot \pmb{n} < 0} \rho(\Theta')~ \phi_s^{\rm ex}(\pmb{r},\pmb{\Omega'},\omega) ~\Big{|} \pmb{\Omega'} \cdot \pmb{n} \Big{|} ~{\rm d}{\Omega'} ~ {\rm{for}}~ (\pmb{r},\pmb\Omega) \in \Gamma^+, \nonumber \end{align} \tag{ 38 }$$

with cos $\Theta'= -\pmb{\Omega'} \cdot \pmb{n}$ and cos $\Theta= \pmb{\Omega} \cdot \pmb{n}$.

The adjoint boundary conditions in the case of specular reflection are:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{BC_adjoint_equations_phi_m_specular_reflection} \phi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)= \rho(\Theta_{\rm sp})~ \phi^{\rm em}(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega)~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma^+ \setminus \Gamma_c^+, \nonumber \end{align} \tag{ 39 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{BC_adjoint_equations_phi_m_specular_reflection_bis} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)= -\big [1-\rho(\Theta)\big ]~(R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)+\rho(\Theta_{\rm sp})~ \phi^{\rm em}(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega)~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma_c^+, \nonumber \end{align} \tag{ 40 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{BC_adjoint_equations_phi_s_specular_reflection} \phi_s^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)= \rho(\Theta_{\rm sp})~ \phi_s^{\rm ex}(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega)~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma^+, \nonumber \end{align} \tag{ 41 }$$

with cos $\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}$ and cos $\Theta= \pmb{\Omega} \cdot \pmb{n}$.

The boundary conditions for the adjoint excitation collimated radiance (related to the third equation of system (33)) fulfill the condition:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{BC_adjoint_equations_phi_c} \phi_c^{\rm ex}(\pmb{r},\omega)= 0 ~ {\rm{for}} ~ \pmb{r} \in \partial \mathcal{D}. \nonumber \end{align} \tag{ 42 }$$

It should be noted that the adjoint variables $\widetilde{\phi_c^{\rm ex}}, \widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}}$ related to the boundary conditions are not used in the adjoint equations model, since these adjoint variables vanish (see appendix D). Moreover, it can be seen that the adjoint equations model takes a similar form to the forward model. The adjoint equations model can be solved in a similar manner to that used to solve the forward model. It also shows that the set of equations should be solved in the order presented in (33). The adjoint emission equation must be solved before the adjoint excitation equation. Thus, the emission adjoint field $\phi^{\rm em}$ obtained first, is used as the source term in the other two equations. After the first equation is solved, the second can be worked out and the excitation diffuse adjoint field $\phi_s^{\rm ex}$ deduced. The resulting $\phi_s^{\rm ex}$ field is used as the source term for the excitation collimated adjoint field $\phi_c^{\rm ex}$.

3.4. Gradient of the objective function

As the adjoint variables $\widetilde{\phi_c^{\rm ex}}, \widetilde{\phi_s^{\rm ex}}, \widetilde{\phi^{\rm em}}$ vanish from the previous remark, the Lagrangian can be written without these variables. The differentiation of the Lagrangian with respect to θ in direction $\delta \theta$ satisfies:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \frac{\partial \mathcal{L}(\theta,\psi_c^{\rm ex},\psi_s^{\rm ex},\psi^{\rm em},\phi_c^{\rm ex},\phi_s^{\rm ex},\phi^{\rm em})}{\partial \theta} \delta \theta=\frac{\partial J(\theta)}{\partial \theta} \delta \theta\nonumber \\ +\Big \langle\phi_c^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_c^{\rm ex}}{\partial \theta}\delta \theta \Big \rangle_S +\Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \theta}\delta \theta \Big \rangle_{SA} +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}} {\partial \theta} \delta \theta \Big \rangle_{SA}.\label{differentiation_Lagrangian} \nonumber \end{align} \tag{ 43 }$$

As the objective function does not depend explicitly on θ (see equation (19)), then we have $\frac{\partial J(\theta)}{\partial \theta}=0$. Using equation (27), equation (43) is reduced to:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Gradient_expressions_1} \Big \langle\pmb\nabla J (\theta) \Big {|} \delta \theta \Big \rangle_S= \Big \langle\phi_c^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_c^{\rm ex}}{\partial \theta}\delta \theta \Big \rangle_S +\Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \theta}\delta \theta \Big \rangle_{SA} +\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}} {\partial \theta} \delta \theta \Big \rangle_{SA}. \nonumber \end{align} \tag{ 44 }$$

This is the expression that evaluates the gradient of the objective function. Applying equation (44) to $\theta=\mu_a^{\rm ex \to em}$ we deduce that:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle\pmb\nabla J (\mu_a^{\rm ex \to em}) \Big {|} \delta \mu_a^{\rm ex \to em} \Big \rangle_S =\Big \langle\phi_c^{\rm ex} \Big{|} \psi_c^{\rm ex} \delta \mu_a^{\rm ex \to em} \Big \rangle_S +\Big \langle\phi_s^{\rm ex} \Big{|} \psi_s^{\rm ex} \delta \mu_a^{\rm ex \to em} \Big \rangle_{SA}\nonumber \\ ~~~~~~~~~~~~~~~~~~~-\Big \langle\phi^{\rm em} \Big{|} \Big (\frac{1}{2\pi}~\frac{\eta}{1-{\rm i}~\omega~\tau}~\Phi^{\rm ex} \Big) \delta \mu_a^{\rm ex \to em} \Big \rangle_{SA}.\label{Gradient_expressions_3} \nonumber \end{align} \tag{ 45 }$$

Then, the components of the objective function gradient, with respect to $\mu_a^{\rm ex \to em}$, are obtained with:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \newcommand{\R}{{\bf R}} \displaystyle \label{Gradient_expressions_4} \pmb\nabla J (\mu_a^{\rm ex \to em}) =\Re ~\Big \{\phi_c^{\rm ex}~\overline{\psi_c^{\rm ex}} + \Big \langle \phi_s^{\rm ex} \Big{|} ~\overline{\psi_s^{\rm ex}} \Big \rangle_A - \frac{1}{2\pi}\Big \langle \phi^{\rm em} \Big{|} \frac{\eta}{1+{\rm i}~\omega~\tau}~\overline{\Phi^{\rm ex}}\Big \rangle_A \Big \} \nonumber \end{align} \tag{ 46 }$$

where $\newcommand{\R}{{\bf R}} \Re ~\{z\}$ is the real part of the complex number z. The right-hand side of (45) is a complex value but only its real part is taken for the gradient of the objective function to be a real value.

3.5. Implementation of the reconstruction algorithm

We present numerical results obtained with the reconstruction algorithm developed in the theoretical section. The calculations were carried out with an Intel Xeon Processor E5-2643, 3.3GHz, 32GoRAM, 8 cores. This last uses Hyper-Threading and Intel C compiler. The computation time for a reconstruction is about 3 hours. The algorithm was used to localize and quantify one and two fluorescent inclusions embedded in a biological medium. The simulations were carried out in a 2D circular domain with a 2 cm radius. The center of the circle is defined as the origin. The chosen size mimics fluorescence tomographic problems that are typical for small animal imaging. The same unstructured triangular mesh composed of 2129 nodal points and 4128 triangular elements was used for the forward and inverse solutions. Each inclusion is a disk of 0.4 cm radius and its center is located at different positions. The intrinsic optical properties of the medium, namely $\mu_a^{\rm ex} = \mu_a^{\rm em} = 0.1~ {\rm{cm}}^{-1}$ and $\mu_s^{\rm ex} = \mu_s^{\rm em} = 100~ {\rm{cm}}^{-1}$ were identical at the excitation and emission wavelengths. The anisotropy factor of the H–G phase function g was set to 0.9. This corresponds to a highly forward-scattering medium. The optical properties are typical for biological tissues in the near-infrared spectral range. The refractive indices of the medium and the surrounding (air) were chosen as n  =  1.4 and $n_{{\rm{out}}} = 1.0$, respectively. The medium surface was assumed to be semi-transparent with specular Fresnel reflection at the interface. ICG (indocyanine green) [10, 37, 84, 85] is the most widely-used fluorophore probe for biomedical applications and was effectively used here. The quantum yield and lifetime of the probe were homogeneously distributed with $\newcommand{\e}{{\rm e}} \eta = 0.012$ and $\tau = 0.52~ {\rm{ns}}$, respectively. We intended to reconstruct the fluorophore absorption coefficient. To this end we allowed determination of the relative ICG concentration in the medium. This process is known as the absorption-contrast mode in FDOT. The fluorescent inclusions represent heterogeneities in the $\mu_a^{\rm ex \to em}$ coefficient [10]. The medium is not ideal, meaning that some ICG molecules remained in the homogeneous background medium leading to residual fluorescence. Therefore, the fluorescence signal detected on the surface had not only been caused by fluorescent inclusions. The fluorophore absorption coefficient in the background was assigned the value $\mu_a^{\rm ex \to em} = 0.01~ {\rm{cm}}^{-1}$. This value was used as the initial estimate in the reconstruction algorithm for all cases described below. A collimated source using ten boundary nodal points (positions all around the disc) was used. These were equally spaced by about 0.1 cm and their centered points were denoted $\pmb{P_c}$. 64 detectors were used that where equally spaced along the boundary nodal points of reflectance. The intensity of the source was modulated at 100 MHz. The use of this modulation frequency leads to a large observation time which is longer than the fluorescence life time of ICG in order to obtain sufficient fluorescence decay. To solve the RTE, the angular space (2πSr) was uniformly subdivided into 32 control solid angles. Furthermore, there were 8 subdivisions into a control solid angle for phase function normalization [76]. The reconstruction was stopped when the normalized difference between two subsequent error functions was smaller than 10⁻⁵. We also tested other values for the initial estimate, namely $\mu_a^{\rm ex \to em} = 10^{-5}~ {\rm{cm}}^{-1}$, $\mu_a^{\rm ex \to em} = 0.03~ {\rm{cm}}^{-1}$ and we obtained reconstruction results close to those for $\mu_a^{\rm ex \to em} = 0.01~ {\rm{cm}}^{-1}$. In the first case, a fluorescent inclusion was embedded inside the homogeneous medium. Its center was located at the exact position (1, 1) as indicated by the dashed circle in figure 1. The true fluorophore absorption coefficient in the inclusion was assigned $\mu_a^{\rm ex \to em} = 0.05~ {\rm{cm}}^{-1}$. For all reconstructions, the synthetic fluorescence data were generated by running the forward model. The true spatial distribution of the fluorophore absorption coefficient we wanted to reconstruct in such a case was used. Particular attention was paid to the optimal configuration of source positions and their collimated directions. Figures 1 and 2(a) depict the results when the collimated directions of the sources were $\pmb{\Omega_c^1}=(0, -1)$ (with $\pmb{P_c^1}=$ (0, 2)) and $\pmb{\Omega_c^2}=({\rm{cos}} (\frac {5\pi}{4}), {\rm{sin}} (\frac {5\pi}{4}))$ (with $\pmb{P_c^2} = (2~{\rm{cos}} (\frac {\pi}{4}), 2~{\rm{sin}} (\frac {\pi}{4}))$). The algorithm gave a well-defined spatial position of the fluorescent inclusion for both directions. However, the $(\pmb{\Omega_c^2}, \pmb{P_c^2})$ source led to a superior reconstruction quality in terms of estimation and localization of the inclusion. With the $(\pmb{\Omega_c^1}, \pmb{P_c^1})$ source, the inclusion was attached to the medium surface and was rather bean shaped. Absorption by the retrieved fluorophore in the inclusion was slightly underestimated relative to its real value. The surface near the inclusion could be illuminated with $(\pmb{\Omega_c^2}, \pmb{P_c^2})$. Then, the circular shape and spatial position were clearly improved as was accuracy. Moreover, the local value retrieved in the inclusion was accurately estimated $(\mu_a^{\rm ex \to em} = 0.05~ {\rm{cm}}^{-1})$. We could then deduce the best configuration to accurately localize and quantify a fluorescent inclusion. It was to illuminate the medium surface near the inclusion using a centered collimated source. The reconstructed images of figures 1 and 2(a) were obtained after 129 and 138 iterations, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the second case we assessed the effect of the inclusion depth on the reconstruction. For this the center of the inclusion can take different exact positions, as shown by the dashed circle in figure 2. The $\pmb{\Omega_c}=({\rm{cos}} (\frac {5\pi}{4}), {\rm{sin}} (\frac {5\pi}{4}))$ direction (with $\pmb{P_c} = (2~{\rm{cos}} (\frac {\pi}{4}), 2~{\rm{sin}} (\frac {\pi}{4}))$) was used as a collimated direction of the source. The estimation and spatial localization of the fluorescent inclusion became less accurate as the depth increased (figure 2). The algorithm failed to reconstruct deeper inclusions, leading to poor image quality. It should be noted that illuminating with more than one collimated sources did not lead to improve the reconstruction. This could be explained by the ill-posedness of the whole problem, in particular light attenuation in the medium. The reconstructed image quality is assessed through the relative RMSE (root mean square error). This is defined between the true and reconstructed values of the fluorophore absorption coefficient:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{error} E~ =~100 \cdot \frac{\big{|} \big{|} \mu_a^{\rm ex \to em}({\rm{reconstructed}})-\mu_a^{\rm ex \to em}({\rm{true}}) \big{|} \big{|}_{2} }{\big{|} \big{|} \mu_a^{\rm ex \to em}({\rm{true}}) \big{|} \big{|}_{2} }, \nonumber \end{align} \tag{ 47 }$$

where $\newcommand{\bi}{\boldsymbol}\big{|} \big{|} \cdot \big{|} \big{|}_2$ denotes Euclidian norm. The relative RMSE was computed over the target region and over the whole reconstructed image domain. As expected, table 1 shows that the relative RMSE increases with the inclusion depth, much more for the target region than that for the whole reconstructed image domain.

As a third case, we assessed the effect of both size and depth of the inclusion on the reconstruction. The relative RMSE over the target region is given in figure 3(a). For each inclusion size, the error decreases as the center of the inclusion moves toward the sources positions. This is because the reconstruction (of shapes and values) is improved when the inclusion is located nearer to the illuminated medium boundary. At each position of the inclusion in the medium, the relative RMSE increases when the inclusion size decreases. As expected, the algorithm failed more to reconstruct a small inclusion located deeper in the medium. Figure 3(b) shows the relative RMSE over the whole reconstructed image domain versus the depth and size of the inclusion. The error is lower (less than 36.5%) when the inclusion size is smallest (R  =  0.2 cm) whatever its position in the medium. When the inclusion is deeper in the medium (located at the center of the medium) the relative RMSE is inverted compared to that over the target region. This means that the reconstruction quality of the whole domain is worst when the inclusion size is bigger and the inclusion is deeper in the medium. Moreover, when R is bigger than 0.2 cm, the relative RMSE over the whole domain increases as the inclusion goes deeper. This is mainly due to the wrong inclusion location inside the homogeneous medium which leads, therefore, to a poor reconstruction quality of the whole domain.

Zoom In Zoom Out Reset image size

In the thourth case, a second inclusion with a lower fluorophore absorption coefficient ($\mu_a^{\rm ex \to em} = 0.03~ {\rm{cm}}^{-1}$) was considered. This allowed us to distinguish between two fluorescent inclusions at different separations. For this purpose, two examples are presented. Figure 4(a) shows the reconstructed image when both inclusions were closed. The centers of the inclusions were located at ($-1, 1$) and ($1, 0.8$) and the phantom medium was simultaneously illuminated by $\pmb{\Omega_c^l}=({\rm{cos}} (\frac {7\pi}{4}), {\rm{sin}} (\frac {7\pi}{4}))$ with $\pmb{P_c^l} = (2~{\rm{cos}} (\frac {3\pi}{4}), 2~{\rm{sin}} (\frac {3\pi}{4}))$ and $\pmb{\Omega_c^r}=({\rm{cos}} (\frac {5\pi}{4}), {\rm{sin}} (\frac {5\pi}{4}))$ with $\pmb{P_c^r} = (2~{\rm{cos}} (\frac {\pi}{4}-\frac {\pi}{32}), 2~{\rm{sin}} (\frac {\pi}{4}-\frac {\pi}{32}))$. Figure 4(b) depicts the result for a somewhat higher separation between the two inclusions. Here, the centers of the inclusions were located at (−1, 1) and (1.4, 0) and the medium was illuminated by $\pmb{\Omega_c^l}=({\rm{cos}} (\frac {7\pi}{4}), {\rm{sin}} (\frac {7\pi}{4}))$ with $\pmb{P_c^l} = (2~{\rm{cos}} (\frac {3\pi}{4}), 2~{\rm{sin}} (\frac {3\pi}{4}))$ and $\pmb{\Omega_c^r}=(-1, 0)$ with $\pmb{P_c^r} = (2, 0)$. The dashed circles in the figures indicate the exact positions of the inclusions. For both separations, the inclusions were accurately recovered and localized in the medium. The algorithm can detect and reconstruct separately the two different fluorescent inclusions. This is so even when they are only separated by small distances (see figure 4(a)). The mean retrieved values of both inclusions in the small and large separation cases are (0.024$~{\rm{cm}}^{-1}$, 0.031$~{\rm{cm}}^{-1}$) and (0.025$~{\rm{cm}}^{-1}$, 0.038$~{\rm{cm}}^{-1}$), respectively. The mean retrieved values of inclusions were always underestimated compared to their exact values. Also, when the inclusions are relatively close to each other the fluorescent photons interfere.

Zoom In Zoom Out Reset image size

In the last case study, the robustness of the reconstruction algorithm was tested in the presence of noisy data with Gaussian distribution [86]. The noise was added to the complex simulated measurement reflectance (module and phase shift). The simulated measurement data are corrupted by adding random errors: 3%, 6% and 10%. Like the first case, the probed medium contains one inclusion. The relative RMSE versus the noise level is shown in table 2. The example of 10% noise is illustrated in figure 5. This figure shows that despite the relatively high noise level, the inclusion is accurately reconstructed. Furthermore, the local retrieved value is well-estimated. This confirmed that the proposed algorithm is robust and efficient, even in the presence of noisy data. However, artifacts and some perturbations of the medium boundary are note compared to the noise-less reconstructed image (see figure 2(a)). This is because detection occurred on the surface and the algorithm converged around a random value. Higher noise levels for the boundary data led to degraded image quality. Relative RMSE decreases with the noise level, much more for the whole reconstructed image domain than that for the target region (see table 2).

Zoom In Zoom Out Reset image size

Table 2. The relative RMSE versus the noise level.

Noise	$\sigma_m$  =  0 (%)	$\sigma_m$  =  3 (%)	$\sigma_m$  =  6 (%)	$\sigma_m = 10$ (%)
Relative RMSE computed over the target region	27.22	29.69	30.97	31.16
Relative RMSE computed over the whole reconstructed image domain	36.14	39.58	49.31	63.35

We studied the inverse problem of fluorescent source. This is based on ICG for a 2D tissue-like phantom subjected to a collimated laser beam. For the first time, we derived the adjoint model for fluorescence molecular imaging based on the RTE. By using the adjoint model an explicit expression of the objective function gradient is obtained, which can be computed efficiently. This is ensured by solving an additional (adjoint) equation for the adjoint variable. The latters's computational cost is equivalent to that of the forward model. The coupled forward and adjoint models were solved with the same modified finite volume method. This had been demonstrated in a previous publication to be highly accurate [76]. The objective function was iteratively minimized using a nonlinear optimization L-BFGS algorithm. Reconstructed images were obtained. These had simulated data for fluorescence signals on the tissue surface. Thus, the spatial fluorophore absorption distribution was assessed taking into account the residual fluorescence in the medium. Positions and directions of the sources were measured on a single fluorescent inclusion. We showed that illuminating the tissue surface from a collimated centered direction near the inclusion gave a better reconstruction quality. Furthermore, we analyzed the effect of both size and depth of the inclusion on the reconstruction. We found that the algorithm failed to reconstruct smaller or deeper inclusions. This was due to light attenuation in the medium. Also, two closely positioned inclusions could be accurately localized. Additionally, their fluorophore absorption coefficients could be quantified. Reconstructions with noisy data were achieved with a reasonable accuracy for several random noise levels. The sum of our results demonstrated that the algorithm is robust. Also, it yields promising results in fluorescence molecular imaging. The present work was a prerequisite study to evaluate the potential of the algorithm. We plan to extend it to 3D geometries (for real applications) using parallel computing with MPI and Open MP, running on a set of multi-core machines.

From (8), (10) and (11), we have:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{calcul_8} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}}=-\frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1-{\rm i}~\omega~\tau(\pmb{r})}. \nonumber \end{align} \tag{ A.1 }$$

The equations (17) and (A.1) yield

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{calcul_9} \Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{SA} =-\int_{\mathcal{D}} \int_{\Omega=2\pi} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)} ~\frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1-{\rm i}~\omega~\tau(\pmb{r})} ~\delta \psi_s^{\rm ex}(\pmb{r},\omega)~{\rm d}\Omega~{\rm d}r, \nonumber \end{align} \tag{ A.2 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{calcul_10} \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_s^{\rm ex}}\Big)^* \phi^{\rm em} =-\frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1+{\rm i}~\omega~\tau(\pmb{r})}~\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega). \nonumber \end{align} \tag{ A.3 }$$

From (4) and (5), we have:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{calcul_2} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}}=-\mu_s^{\rm ex}(\pmb{r})~p^{\rm ex}(\pmb{\Omega_c}\cdot\pmb{\Omega}). \nonumber \end{align} \tag{ A.4 }$$

The equations (17) and (A.4) yield

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \Big \langle\phi_s^{\rm ex} \Big{|} \frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA} = \int_{\mathcal{D}} -\mu_s^{\rm ex}(\pmb{r})~\Big (\int_{\Omega=2\pi} p^{\rm ex}(\pmb{\Omega_c}\cdot\pmb{\Omega})~\overline{\phi_s^{\rm ex}(\pmb{r},\pmb\Omega,\omega)}~{\rm d}\Omega \Big) ~\delta \psi_c^{\rm ex}(\pmb{r},\omega)~{\rm d}r, \label{calcul_3} \nonumber \end{align} \tag{ A.5 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{calcul_4} \Big(\frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}}\Big)^{\overline{*}} \phi_s^{\rm ex}=-\mu_s^{\rm ex}(\pmb{r})~ \int_{\Omega=2\pi} p^{\rm ex}(\pmb{\Omega_c}\cdot\pmb{\Omega}) ~\phi_s^{\rm ex}(\pmb{r},\pmb\Omega,\omega)~{\rm d}\Omega ~ {\rm{with}} ~ \frac{1}{2 \pi} \Big(\frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}}\Big)^{\overline{*}} =\Big(\frac {\partial \mathcal{R}_s^{\rm ex}}{\partial \psi_c^{\rm ex}}\Big)^*. \nonumber \end{align} \tag{ A.6 }$$

From (8), (10) and (11), we have:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{calcul_5} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}}=-\frac{1}{2\pi}~\frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1-{\rm i}~\omega~\tau(\pmb{r})}. \nonumber \end{align} \tag{ A.7 }$$

The equations (17) and (A.7) yield

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber &\Big \langle\phi^{\rm em} \Big{|} \frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_{SA} =-\frac{1}{2\pi}~\int_{\mathcal{D}} \frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1-{\rm i}~\omega~\tau(\pmb{r})} \Big (~ \int_{\Omega=2\pi} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)} ~{\rm d}\Omega \Big) ~\delta \psi_c^{\rm ex}(\pmb{r},\omega)~{\rm d}r\nonumber \\ &\qquad \qquad \qquad \qquad =\Big \langle \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}}\Big)^{\overline{*}} \phi^{\rm em} \Big{|} ~\delta \psi_c^{\rm ex} \Big \rangle_S, \label{calcul_6} \nonumber \end{align} \tag{ A.8 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{calcul_7} \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}}\Big)^{\overline{*}} \phi^{\rm em} =-\frac{1}{2\pi}~\frac{\eta~\mu_a^{\rm ex \to em}(\pmb{r})}{1+{\rm i}~\omega~\tau(\pmb{r})} \int_{\Omega=2\pi} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~{\rm d}\Omega~ {\rm{and}} ~ \frac{1}{2 \pi} \Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}}\Big)^{\overline{*}} =\Big(\frac {\partial \mathcal{R}^{\rm em}}{\partial \psi_c^{\rm ex}}\Big)^*. \nonumber \end{align} \tag{ A.9 }$$

The boundary conditions of the adjoint RTE in the time domain for semi-transparent boundaries (with diffuse reflection) were derived by Boulanger et al [69]. Here we propose to adapt these boundary conditions for the adjoint RTE in the frequency domain. Let ψ and ϕ be two functions of the space $L^2(\mathcal{D} \times \Sigma)$. Using integration by parts and the divergence theorem, we obtain:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Adjoint_transport_operator_1} \Big \langle \phi \Big{|} \big(\pmb\Omega \cdot \pmb\nabla \big) \psi \Big \rangle_{SA} = \int_{\partial \mathcal{D}} \int_{\Omega=2\pi} \overline{\phi(\pmb{r},\pmb\Omega,\omega)}~\psi(\pmb{r},\pmb\Omega,\omega)~\pmb\Omega \cdot \pmb n ~{\rm d}\Omega~{\rm d}r -\Big\langle\big(\pmb\Omega \cdot \pmb\nabla \big) \phi \Big{|} \psi \Big \rangle_{SA}. \nonumber \end{align} \tag{ B.1 }$$

Then $\newcommand{\bi}{\boldsymbol}-\big(\pmb\Omega \cdot \pmb\nabla \big)$ is the adjoint of the transport operator $\newcommand{\bi}{\boldsymbol}\big(\pmb\Omega \cdot \pmb\nabla \big)$, provided the boundary term vanishes for all $\psi \in L^2(\mathcal{D} \times \Sigma)$, i.e.:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_2} \int_{\partial \mathcal{D}} \int_{\pmb\Omega' \cdot \pmb n <0} \overline{\phi(\pmb{r},\pmb\Omega',\omega)}~\psi(\pmb{r},\pmb\Omega',\omega) ~\pmb\Omega' \cdot \pmb n~{\rm d}\Omega'~{\rm d}r =-\int_{\partial \mathcal{D}} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi(\pmb{r},\pmb\Omega,\omega)}~\psi(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r. \nonumber \end{align} \tag{ B.2 }$$

If $\partial \mathcal{D}$ is assumed to be transparent (no reflection at the boundary), the left-hand side of (B.2) is equal to zero since $\psi(\pmb{r}, \pmb\Omega, \omega)=0$ for $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} (\pmb{r}, \pmb\Omega) \in \Gamma^-$. This implies that the right-hand side of (B.2) is also equal to zero and ϕ must satisfy the adjoint boundary conditions $\phi(\pmb{r}, \pmb\Omega, \omega)=0$ for $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} (\pmb{r}, \pmb\Omega) \in \Gamma^+$. If $\partial \mathcal{D}$ is assumed to be semi-transparent, using the boundary conditions (7) in the case of diffuse reflection, the left-hand side of (B.2) is equal to:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Adjoint_transport_operator_3} \int_{\partial \mathcal{D}} \Big(\int_{\pmb\Omega' \cdot \pmb n <0} \overline{\phi(\pmb{r},\pmb\Omega',\omega)}~\pmb\Omega' \cdot \pmb n~{\rm d}\Omega' \Big)~\Big(\frac{1}{\pi}\int_{\pmb{\Omega} \cdot \pmb{n} > 0} \rho(\Theta)~ \psi(\pmb{r},\pmb{\Omega},\omega)~ \pmb{\Omega} \cdot \pmb{n} ~{\rm d}{\Omega} \Big) ~{\rm d}r ~ {\rm{with}}~ {\rm{cos}}~\Theta= \pmb{\Omega} \cdot \pmb{n}, \nonumber \end{align} \tag{ B.3 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Adjoint_transport_operator_3_discrete_sum} =\frac{1}{\pi}\int_{\partial \mathcal{D}} \Big(\underset{N_l \to \infty}{{\rm{lim}}} \sum_{\big\{l,~\pmb{\Omega^l} \cdot \pmb{n} < 0 \big\}}^{N_l} \overline{\phi(\pmb{r},\pmb\Omega^l,\omega)}~ \pmb{\Omega^l} \cdot \pmb{n} ~\omega^l \Big) \Big(\underset{N_k \to \infty}{{\rm{lim}}} \sum_{\big\{k,~\pmb{\Omega^k} \cdot \pmb{n} > 0 \big\}}^{N_k} \rho(\Theta^k)~ \psi(\pmb{r},\pmb{\Omega^k},\omega)~ \pmb{\Omega^k} \cdot \pmb{n} ~\omega^k \Big) ~{\rm d}r, \nonumber \end{align} \tag{ B.4 }$$

where N_l and N_k are the total numbers of discrete ordinates $\pmb\Omega^l$ and $\pmb\Omega^k$, respectively. The quadrature weights w_l and w_k are associated respectively to the discrete ordinates $\pmb\Omega^l$ and $\pmb\Omega^k$ [10]. The discrete angles and directions are chosen such that:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{relation_angles} \Theta^k=\Theta^l \in \Sigma, ~ {\rm{cos}}~ \Theta^k= \pmb{\Omega^k} \cdot \pmb{n} = \pmb{\Omega^l} \cdot (-\pmb{n}) = {\rm{cos}} ~\Theta^l ~ {\rm{for~all}}~ 1 \leqslant k \leqslant {N_k}~ {\rm{and}} ~ 1 \leqslant l \leqslant {N_l}. \nonumber \end{align} \tag{ B.5 }$$

Then, the relation (B.4) is also equal to:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_3_discrete_sum_bis} \frac{1}{\pi} \int_{\partial \mathcal{D}} \underset{N_k \to \infty}{{\rm{lim}}} \underset{N_l \to \infty}{{\rm{lim}}} \sum_{\big\{k,~\pmb{\Omega^k} \cdot \pmb{n} > 0 \big\}}^{N_k} \sum_{\big\{l,~\pmb{\Omega^l} \cdot \pmb{n} < 0 \big\}}^{N_l} \rho(\Theta^k)~\overline{\phi(\pmb{r},\pmb\Omega^l,\omega)}~ \psi(\pmb{r},\pmb{\Omega^k},\omega)~ \pmb{\Omega^k} \cdot \pmb{n}~\pmb{\Omega^l} \cdot \pmb{n} ~\omega^k~\omega^l ~{\rm d}r, \nonumber \end{align} \tag{ B.6 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_3_discrete_sum_bis_bis} =\frac{1}{\pi}\int_{\partial \mathcal{D}} \underset{N_k \to \infty}{{\rm{lim}}} \underset{N_l \to \infty}{{\rm{lim}}} \sum_{\big\{k,~\pmb{\Omega^k} \cdot \pmb{n} > 0 \big\}}^{N_k} \sum_{\big\{l,~\pmb{\Omega^l} \cdot \pmb{n} < 0 \big\}}^{N_l} \rho(\Theta^l)~\overline{\phi(\pmb{r},\pmb\Omega^l,\omega)}~ \psi(\pmb{r},\pmb{\Omega^k},\omega)~ \pmb{\Omega^k} \cdot \pmb{n}~\pmb{\Omega^l} \cdot \pmb{n} ~\omega^k~\omega^l ~{\rm d}r, \nonumber \end{align} \tag{ B.7 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Adjoint_transport_operator_3_discrete_sum_bis_bis_bis} =\frac{1}{\pi} \int_{\partial \mathcal{D}} \underset{N_k \to \infty}{{\rm{lim}}} \sum_{\big\{k,~\pmb{\Omega^k} \cdot \pmb{n} > 0 \big\}}^{N_k} \Big(\underset{N_l \to \infty}{{\rm{lim}}} \sum_{\big\{l,~\pmb{\Omega^l} \cdot \pmb{n} < 0 \big\}}^{N_l} \rho(\Theta^l)~\overline{\phi(\pmb{r},\pmb\Omega^l,\omega)}~ \pmb{\Omega^l} \cdot \pmb{n} ~\omega^l \Big) \psi(\pmb{r},\pmb{\Omega^k},\omega)~ \pmb{\Omega^k} \cdot \pmb{n} ~\omega^k ~{\rm d}r, \nonumber \end{align} \tag{ B.8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Adjoint_transport_operator_4} = \int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \Big(\frac{1}{\pi}\int_{\pmb\Omega' \cdot \pmb n <0} \rho(\Theta')~ \overline{\phi(\pmb{r},\pmb\Omega',\omega)}~\pmb\Omega' \cdot \pmb n~{\rm d}\Omega' \Big)~ \psi(\pmb{r},\pmb{\Omega},\omega)~ \pmb{\Omega} \cdot \pmb{n} ~{\rm d}{\Omega}~{\rm d}r ~ {\rm{with}}~ {\rm{cos}}~\Theta'= -\pmb{\Omega'} \cdot \pmb{n}. \nonumber \end{align} \tag{ B.9 }$$

By identification of the right-hand side of (B.2) with (B.9), the relation yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Adjoint_transport_operator_5} \phi(\pmb{r},\pmb\Omega,\omega)=\frac{1}{\pi}\int_{\pmb\Omega' \cdot \pmb n <0} \rho(\Theta')~\phi(\pmb{r},\pmb\Omega',\omega)~\Big{|}\pmb\Omega' \cdot \pmb n \Big{|}~{\rm d}\Omega', \nonumber \end{align} \tag{ B.10 }$$

for $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} (\pmb{r}, \pmb\Omega) \in \Gamma^+$ with cos $\Theta'= -\pmb{\Omega'} \cdot \pmb{n}$.

Using the boundary conditions (7) in the case of specular reflection, the left-hand side of (B.2) is now equal to:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_6} \int_{\partial \mathcal{D}} \int_{\pmb\Omega' \cdot \pmb n <0} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \overline{\phi(\pmb{r},\pmb\Omega',\omega)}~\pmb\Omega' \cdot \pmb n~ \rho(\Theta_{\rm sp})~ \psi(\pmb{r},\pmb{\Omega},\omega)~ \delta(\pmb{\Omega}-\pmb{\Omega_{\rm sp}}) ~{\rm d}{\Omega}~{\rm d}\Omega' ~{\rm d}r, \nonumber \end{align} \tag{ B.11 }$$

where cos $\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}>0$ with $\pmb{\Omega_{\rm sp}}=\pmb\Omega'-2(\pmb{\Omega'} \cdot \pmb{n})~\pmb{n}$. In this case, $\pmb{\Omega'} \cdot \pmb{n}=-\pmb{\Omega_{\rm sp}} \cdot \pmb{n}$ and (B.11) is equal to:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_7} -\int_{\partial \mathcal{D}} \int_{\pmb\Omega' \cdot \pmb n <0} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \rho(\Theta_{\rm sp})~ \overline{\phi(\pmb{r},\pmb\Omega',\omega)}~\pmb{\Omega} \cdot \pmb n ~ \psi(\pmb{r},\pmb{\Omega},\omega)~ \delta(\pmb{\Omega}-\pmb{\Omega_{\rm sp}}) ~{\rm d}{\Omega}~{\rm d}\Omega' ~{\rm d}r, \nonumber \end{align} \tag{ B.12 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_7_bis} =-\int_{\partial \mathcal{D}} \underset{N_k \to \infty}{{\rm{lim}}} \underset{N_l \to \infty}{{\rm{lim}}} \sum_{\big\{k,~\pmb{\Omega^k} \cdot \pmb{n} > 0 \big\}}^{N_k} \sum_{\big\{l,~\pmb{\Omega^l} \cdot \pmb{n} < 0 \big\}}^{N_l} \rho(\Theta_{\rm sp})~\overline{\phi(\pmb{r},\pmb\Omega^l,\omega)}~\pmb{\Omega^k} \cdot \pmb n ~ \psi(\pmb{r},\pmb{\Omega^k},\omega)~ \delta(\pmb{\Omega^k}-\pmb{\Omega_{\rm sp}}) ~{\rm d}r. \nonumber \end{align} \tag{ B.13 }$$

Using (B.5), it follows that:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_7_bis_bis} =-\int_{\partial \mathcal{D}} \underset{N_k \to \infty}{{\rm{lim}}} \underset{N_l \to \infty}{{\rm{lim}}} \sum_{\big\{k,~\pmb{\Omega^k} \cdot \pmb{n} > 0 \big\}}^{N_k} \sum_{\big\{l,~\pmb{\Omega^l} \cdot \pmb{n} < 0 \big\}}^{N_l} \rho(\Theta_{\rm sp})~\overline{\phi(\pmb{r},\pmb\Omega^l,\omega)}~\pmb{\Omega^k} \cdot \pmb n ~ \psi(\pmb{r},\pmb{\Omega^k},\omega)~ \delta(\pmb{\Omega^l}+\pmb{\Omega_{\rm sp}}) ~{\rm d}r, \nonumber \end{align} \tag{ B.14 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_8} =-\int_{\partial \mathcal{D}} \int_{\pmb\Omega' \cdot \pmb n <0} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \rho(\Theta_{\rm sp})~ \overline{\phi(\pmb{r},\pmb\Omega',\omega)}~\pmb{\Omega} \cdot \pmb n ~ \psi(\pmb{r},\pmb{\Omega},\omega)~ \delta(\pmb{\Omega'}+\pmb{\Omega_{\rm sp}}) ~{\rm d}{\Omega}~{\rm d}\Omega' ~{\rm d}r, \nonumber \end{align} \tag{ B.15 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_9} =-\int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \rho(\Theta_{\rm sp})~ \overline{\phi(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega)}~ \psi(\pmb{r},\pmb{\Omega},\omega)~\pmb{\Omega} \cdot \pmb n~{\rm d}{\Omega}~{\rm d}r. \nonumber \end{align} \tag{ B.16 }$$

By identification of the right-hand side of (B.2) with (B.16), the relation yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{Adjoint_transport_operator_10} \phi(\pmb{r},\pmb{\Omega},\omega)= \rho(\Theta_{\rm sp})~ \phi(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega) ~ {\rm{for}} ~ (\pmb{r},\pmb\Omega) \in \Gamma^+ ~ {\rm{with}} ~ {\rm{cos}}~\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}. \nonumber \end{align} \tag{ B.17 }$$

Now, let ψ and ϕ be two functions of the space $L^2(\mathcal{D})$. Using integration by parts and divergence theorem, we obtain:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Adjoint_transport_operator_11} \Big \langle \phi \Big{|} \big(\pmb\Omega_c \cdot \pmb\nabla \big) \psi \Big \rangle_S = \int_{\partial \mathcal{D}} \overline{\phi(\pmb{r},\omega)}~\psi(\pmb{r},\omega)~\pmb\Omega_c \cdot \pmb n ~{\rm d}r -\Big\langle\big(\pmb\Omega_c \cdot \pmb\nabla \big) \phi \Big{|} \psi \Big \rangle_S. \nonumber \end{align} \tag{ B.18 }$$

Then $\newcommand{\bi}{\boldsymbol}-\big(\pmb\Omega_c \cdot \pmb\nabla \big)$ is the adjoint of the transport operator $\newcommand{\bi}{\boldsymbol}\big(\pmb\Omega_c \cdot \pmb\nabla \big)$, provided the boundary term vanishes for all $\psi \in L^2(\mathcal{D})$. Using the boundary conditions (2), it follows that:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_13} \int_{\partial \mathcal{D}_c} \overline{\phi(\pmb{r},\omega)}~\Upsilon(\pmb{r},\omega) ~\pmb\Omega_c \cdot \pmb n~{\rm d}r=0. \nonumber \end{align} \tag{ B.19 }$$

Hence, we obtain:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Adjoint_transport_operator_14} \phi(\pmb{r},\omega)= 0 ~ {\rm{almost~everywhere}}~ {\rm{for}} ~ \pmb{r} \in \partial \mathcal{D}_c ~ {\rm{and}} ~ \phi(\pmb{r},\omega)= 0 ~ ~ {\rm{elsewhere}}. \nonumber \end{align} \tag{ B.20 }$$

Let $\psi^{\rm em}$ and $\phi^{\rm em}$ be two functions of the space $L^2(\mathcal{D} \times \Sigma)$. Applying integration by parts and the divergence theorem to the transport operator of the first equation of (30), we obtain:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber &\Big \langle \phi^{\rm em} \Big{|} \big(\pmb\Omega \cdot \pmb\nabla \big) \psi^{\rm em} \Big \rangle_{SA}+\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \psi^{\rm em}) \Big \rangle_B = \int_{\partial \mathcal{D}} \int_{\Omega=2\pi} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega)~\pmb\Omega \cdot \pmb n ~{\rm d}\Omega~{\rm d}r\nonumber \\ &\qquad \qquad \qquad \qquad -\Big\langle\big(\pmb\Omega \cdot \pmb\nabla \big) \phi^{\rm em} \Big{|} \psi^{\rm em} \Big \rangle_{SA}+\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \psi^{\rm em}) \Big \rangle_B.\label{Adjoint_transport_operator_appendix_c} \nonumber \end{align} \tag{ C.1 }$$

To ensure that $\newcommand{\bi}{\boldsymbol}-\big(\pmb\Omega \cdot \pmb\nabla \big)$ is the adjoint of the transport operator $\newcommand{\bi}{\boldsymbol}\big(\pmb\Omega \cdot \pmb\nabla \big)$, the boundary terms of the right-hand side of (C.1) must vanish for all $\psi^{\rm em} \in L^2(\mathcal{D} \times \Sigma)$, i.e.:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \int_{\partial \mathcal{D}} \int_{\pmb\Omega' \cdot \pmb n <0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega',\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega',\omega) ~\pmb\Omega' \cdot \pmb n~{\rm d}\Omega'~{\rm d}r\nonumber \\ +\int_{\partial \mathcal{D}} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r +\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \psi^{\rm em}) \Big \rangle_B=0.\label{new_term_BC_2} \nonumber \end{align} \tag{ C.2 }$$

If $\partial \mathcal{D}$ is assumed to be transparent (no reflection at the boundary), the first term of (C.2) is equal to zero since $\psi^{\rm em}(\pmb{r}, \pmb\Omega, \omega)=0$ for $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} (\pmb{r}, \pmb\Omega) \in \Gamma^-$ and:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{new_term_BC_3} \int_{\partial \mathcal{D}} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r +\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \psi^{\rm em}) \Big \rangle_B=0. \nonumber \end{align} \tag{ C.3 }$$

As (C.3) has to be satisfied for all $\psi^{\rm em} \in L^2(\mathcal{D} \times \Sigma)$, $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} \phi^{\rm em}(\pmb{r}, \pmb\Omega, \omega)=0 ~{\rm{for}}~ (\pmb{r}, \pmb\Omega) \in \Gamma^+ \setminus \Gamma_c^+$ and using (20) it follows that:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{new_term_BC_4} &\int_{\partial \mathcal{D}_c} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r \nonumber \\ & =-\int_{\partial \mathcal{D}_c}\int_{\pmb\Omega \cdot \pmb n >0} (\overline{R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega)~\pmb\Omega \cdot \pmb n ~{\rm d}\Omega~{\rm d}r. \nonumber \end{align} \tag{ C.4 }$$

By identification, we have the following adjoint boundary conditions:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{new_term_BC_5} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)=- (R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma_c^+. \nonumber \end{align} \tag{ C.5 }$$

If $\partial \mathcal{D}$ is assumed to be semi-transparent with diffuse reflection, the first term of (C.2) is equal to (B.9) and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \Big(\frac{1}{\pi}\int_{\pmb\Omega' \cdot \pmb n <0} \rho(\Theta')~ \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega',\omega)}~\pmb\Omega' \cdot \pmb n~{\rm d}\Omega' \Big)~ \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)~ \pmb{\Omega} \cdot \pmb{n} ~{\rm d}{\Omega}~{\rm d}r\nonumber \\ +\int_{\partial \mathcal{D}} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r +\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \psi^{\rm em}) \Big \rangle_B=0\label{new_term_BC_6} \nonumber \end{align} \tag{ C.6 }$$

with cos $\Theta'= -\pmb{\Omega'} \cdot \pmb{n}$. By identification, the following relation yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \newcommand{\B}{{\bf B}} \displaystyle \label{new_term_BC_7} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)=\frac{1}{\pi}\int_{\pmb\Omega' \cdot \pmb n <0} \rho(\Theta')~\phi^{\rm em}(\pmb{r},\pmb\Omega',\omega)~\Big{|}\pmb\Omega' \cdot \pmb n \Big{|}~{\rm d}\Omega'~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma^+ \setminus \Gamma_c^+, \nonumber \end{align} \tag{ C.7 }$$

and using (20) it follows that:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber &\int_{\partial \mathcal{D}_c} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \Big(\frac{1}{\pi}\int_{\pmb\Omega' \cdot \pmb n <0} \rho(\Theta')~ \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega',\omega)}~\pmb\Omega' \cdot \pmb n~{\rm d}\Omega' \Big)~ \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)~ \pmb{\Omega} \cdot \pmb{n} ~{\rm d}{\Omega}~{\rm d}r\nonumber \\ \nonumber &\qquad \qquad \qquad ~+\int_{\partial \mathcal{D}_c} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r\nonumber \\ & \qquad \qquad =-\int_{\partial \mathcal{D}_c}\int_{\pmb\Omega \cdot \pmb n >0} \big [1-\rho(\Theta)\big ]~(\overline{R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega)~\pmb\Omega \cdot \pmb n ~{\rm d}\Omega~{\rm d}r,\label{new_term_BC_8} \nonumber \end{align} \tag{ C.8 }$$

with cos $\Theta= \pmb{\Omega} \cdot \pmb{n}$. By identification, the following relation yields:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \newcommand{\B}{{\bf B}} \displaystyle \label{new_term_BC_9} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)=& -\big [1-\rho(\Theta)\big ]~(R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)\nonumber \\ & +\frac{1}{\pi}\int_{\pmb\Omega' \cdot \pmb n <0} \rho(\Theta')~\phi^{\rm em}(\pmb{r},\pmb\Omega',\omega)~\Big{|}\pmb\Omega' \cdot \pmb n \Big{|}~{\rm d}\Omega'~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma_c^+. \nonumber \end{align} \tag{ C.9 }$$

If $\partial \mathcal{D}$ is assumed to be semi-transparent with specular reflection, the first term of (C.2) is equal to (B.16) and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber -\int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \rho(\Theta_{\rm sp})~ \overline{\phi^{\rm em}(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega)}~ \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)~\pmb{\Omega} \cdot \pmb n~{\rm d}{\Omega}~{\rm d}r\nonumber \\ +\int_{\partial \mathcal{D}} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r +\Big \langle R^{\rm em}-d_{\rm obs}^{\rm em} \Big{|} (H \psi^{\rm em}) \Big \rangle_B=0\label{new_term_BC_10} \nonumber \end{align} \tag{ C.10 }$$

with cos $\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}$. By identification, the following relation yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{new_term_BC_11} \phi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)= \rho(\Theta_{\rm sp})~ \phi^{\rm em}(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega) ~ {\rm{for}} ~ (\pmb{r},\pmb\Omega) \in \Gamma^+ \setminus \Gamma_c^+, \nonumber \end{align} \tag{ C.11 }$$

and using (20):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber -\int_{\partial \mathcal{D}_c} \int_{\pmb{\Omega} \cdot \pmb{n} > 0} \rho(\Theta_{\rm sp})~ \overline{\phi^{\rm em}(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega)}~ \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)~\pmb{\Omega} \cdot \pmb n~{\rm d}{\Omega}~{\rm d}r\nonumber \\ +\int_{\partial \mathcal{D}_c} \int_{\pmb\Omega \cdot \pmb n >0} \overline{\phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega) ~\pmb\Omega \cdot \pmb n~{\rm d}\Omega~{\rm d}r \nonumber \end{align} \tag{ C.12 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle =-\int_{\partial \mathcal{D}_c}\int_{\pmb\Omega \cdot \pmb n >0} \big [1-\rho(\Theta)\big ]~(\overline{R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)}~\psi^{\rm em}(\pmb{r},\pmb\Omega,\omega)~\pmb\Omega \cdot \pmb n ~{\rm d}\Omega~{\rm d}r,\label{new_term_BC_12} \nonumber \end{align} \tag{ C.13 }$$

with cos $\Theta= \pmb{\Omega} \cdot \pmb{n}$. By identification, the following relation yields:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\G}{{{\leftrightarrow}{\bf G}}} \displaystyle \label{new_term_BC_13} \phi^{\rm em}(\pmb{r},\pmb\Omega,\omega)= -\big [1-\rho(\Theta)\big ]~(R^{\rm em}-d_{\rm obs}^{\rm em}) (\pmb{r},\omega)+\rho(\Theta_{\rm sp})~ \phi^{\rm em}(\pmb{r},-\pmb{\Omega_{\rm sp}},\omega)~{\rm{for}}~(\pmb{r},\pmb\Omega) \in \Gamma_c^+. \nonumber \end{align} \tag{ C.14 }$$

These terms are determined from (30).

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \label{Appendix_D_1} \Big \langle \widetilde{\phi_c^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_c^{\rm ex}} }{\partial \psi_c^{\rm ex}} \delta \psi_c^{\rm ex} \Big \rangle_B= \int_{\partial \mathcal{D}_c} \overline{\widetilde{\phi_c^{\rm ex}}(\pmb{r},\omega)}~\delta \psi_c^{\rm ex}(\pmb{r},\omega) ~{\rm d}r=0 \nonumber \end{align} \tag{ D.1 }$$

Then,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Appendix_D_2} \widetilde{\phi_c^{\rm ex}}(\pmb{r},\omega)= 0 ~ {\rm{almost~everywhere}} ~{\rm{for}} ~ \pmb{r} \in \partial \mathcal{D}_c. \nonumber \end{align} \tag{ D.2 }$$

In the case of diffuse reflection, from (7) it yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber &\qquad \qquad \Big \langle \widetilde{\phi_s^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_s^{\rm ex}} }{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{BA^-}= \int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi_s^{\rm ex}}(\pmb{r},\pmb{\Omega},\omega)}~\delta \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)~{\rm d}\Omega~ {\rm d}r\nonumber \\ &-\frac{1}{\pi} \int_{\partial \mathcal{D}}\Bigg(\int_{\pmb{\Omega'} \cdot \pmb{n} > 0} \rho(\Theta')~\delta \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega'},\omega) ~\pmb{\Omega'} \cdot \pmb{n} ~{\rm d}{\Omega'}\Bigg) \Bigg(\int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi_s^{\rm ex}}(\pmb{r},\pmb{\Omega},\omega)}~{\rm d}\Omega \Bigg)~ {\rm d}r=0~ {\rm{with}}~ {\rm{cos}}~\Theta'= \pmb{\Omega'} \cdot \pmb{n}.\label{Appendix_D_3_diffuse_reflection} \nonumber \end{align} \tag{ D.3 }$$

In the case of specular reflection, from (7) it yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle \widetilde{\phi_s^{\rm ex}} \Big{|} \frac {\partial \widetilde{\mathcal{R}_s^{\rm ex}} }{\partial \psi_s^{\rm ex}} \delta \psi_s^{\rm ex} \Big \rangle_{BA^-}= \int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi_s^{\rm ex}}(\pmb{r},\pmb{\Omega},\omega)}~\delta \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega},\omega)~{\rm d}\Omega~ {\rm d}r\nonumber \\ -\int_{\partial \mathcal{D}}\Bigg(\rho(\Theta_{\rm sp})~\delta \psi_s^{\rm ex}(\pmb{r},\pmb{\Omega_{\rm sp}},\omega)\Bigg) \Bigg(\int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi_s^{\rm ex}}(\pmb{r},\pmb{\Omega},\omega)}~{\rm d}\Omega \Bigg)~ {\rm d}r=0~ {\rm{with}}~ {\rm{cos}}~\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}.\label{Appendix_D_3_specular_reflection} \nonumber \end{align} \tag{ D.4 }$$

A solution of equations (D.3) and (D.4) is given by $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} \widetilde{\phi_s^{\rm ex}}(\pmb{r}, \pmb{\Omega}, \omega)= 0 ~ {\rm{for}} ~ (\pmb{r}, \pmb\Omega) \in \Gamma^-.$

In the case of diffuse reflection, from (12) it yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber &\Big \langle \widetilde{\phi^{\rm em}} \Big{|} \frac {\partial \widetilde{\mathcal{R}^{\rm em}} }{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{BA^-}= \int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi^{\rm em}}(\pmb{r},\pmb{\Omega},\omega)}\delta \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)~{\rm d}\Omega~ {\rm d}r\nonumber \\ &-\frac{1}{\pi} \int_{\partial \mathcal{D}}\Bigg(\int_{\pmb{\Omega'} \cdot \pmb{n} > 0} \rho(\Theta')~\delta \psi^{\rm em}(\pmb{r},\pmb{\Omega'},\omega) \pmb{\Omega'} \cdot \pmb{n} ~{\rm d}{\Omega'}\Bigg) \Bigg(\int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi^{\rm em}}(\pmb{r},\pmb{\Omega},\omega)}~{\rm d}\Omega \Bigg) {\rm d}r\!=0~ {\rm{with}}~{\rm{cos}}~\Theta'= \pmb{\Omega'} \cdot \pmb{n}.\label{Appendix_D_4_diffuse_reflection} \nonumber \end{align} \tag{ D.5 }$$

In the case of specular reflection, from (12) it yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\B}{{\bf B}} \displaystyle \nonumber \Big \langle \widetilde{\phi^{\rm em}} \Big{|} \frac {\partial \widetilde{\mathcal{R}^{\rm em}} }{\partial \psi^{\rm em}} \delta \psi^{\rm em} \Big \rangle_{BA^-}= \int_{\partial \mathcal{D}} \int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi^{\rm em}}(\pmb{r},\pmb{\Omega},\omega)}~\delta \psi^{\rm em}(\pmb{r},\pmb{\Omega},\omega)~{\rm d}\Omega~ {\rm d}r\nonumber \\ -\int_{\partial \mathcal{D}}\Bigg(\rho(\Theta_{\rm sp})~\delta \psi^{\rm em}(\pmb{r},\pmb{\Omega_{\rm sp}},\omega)\Bigg) \Bigg(\int_{\pmb{\Omega} \cdot \pmb{n} < 0} \overline{\widetilde{\phi^{\rm em}}(\pmb{r},\pmb{\Omega},\omega)}~{\rm d}\Omega \Bigg)~ {\rm d}r=0~ {\rm{with}}~ {\rm{cos}}~\Theta_{\rm sp}= \pmb{\Omega_{\rm sp}} \cdot \pmb{n}.\label{Appendix_D_4_specular_reflection} \nonumber \end{align} \tag{ D.6 }$$

A solution of equations (D.5) and (D.6) is given by $\newcommand{\G}{{{\leftrightarrow}{\bf G}}} \widetilde{\phi^{\rm em}}(\pmb{r}, \pmb{\Omega}, \omega)= 0 ~ {\rm{for}} ~ (\pmb{r}, \pmb\Omega) \in \Gamma^-.$