On Iterative Algorithms for Quantitative Photoacoustic Tomography in the Radiative Transport Regime

TL;DR

This paper develops and compares iterative algorithms for quantitative photoacoustic tomography using the radiative transfer equation, providing efficient methods for reconstructing optical tissue properties with proven convergence.

Contribution

It introduces an improved fixed-point iterative method and applies the Barzilai-Borwein method for simultaneous reconstruction of optical coefficients in QPAT based on RTE.

Findings

The fixed-point method converges reliably for absorption coefficient reconstruction.

The BB method effectively reconstructs both absorption and scattering coefficients.

Simulation results demonstrate the efficiency of the proposed algorithms.

Abstract

In this paper, we describe the numerical reconstruction method for quantitative photoacoustic tomography (QPAT) based on the radiative transfer equation (RTE), which models light propagation more accurately than diffusion approximation (DA). We investigate the reconstruction of absorption coefficient and/or scattering coefficient of biological tissues. Given the scattering coefficient, an improved fixed-point iterative method is proposed to retrieve the absorption coefficient for its cheap computational cost. And we prove the convergence. To retrieve two coefficients simultaneously, Barzilai-Borwein (BB) method is applied. Since the reconstruction of optical coefficients involves the solution of original and adjoint RTEs in the framework of optimization, an efficient solver with high accuracy is improved from~\cite{Gao}. Simulation experiments illustrate that the improved fixed-point iterative method and the BB method are the comparative methods for QPAT in two cases.

Tables

Table 1. Table 1: Relative L 2 superscript 𝐿 2 L^{2} errors on reconstructions by BB method.

	noise-free data		noisy data
	${ϵ}_{{\mu }_a}$	${ϵ}_{{\mu }_s}$	${ϵ}_{{\mu }_a}$	${ϵ}_{{\mu }_s}$
1	4.61e-2	1.52e-1	1.09e-1	1.81e-1
2	2.65e-2	1.53e-1	7.71e-2	1.64e-1
3	3.17e-2	1.09e-1	1.29e-1	1.37e-1
4	5.32e-2	1.26e-1	1.28e-1	1.49e-1

Equations

$${\Gamma }_{-}:=$$

$${\Gamma }_{+}:=$$

$$\left[\theta \cdot \nabla +{\mu }_a(x)+{\mu }_s(x)\right]\varphi (x,\theta )-{\mu }_s(x){\oint }_{{\mathcal{S}}^{n-1}}f(\theta ,{\theta }^{\prime })\varphi (x,{\theta }^{\prime })\mathrm{d}{\theta }^{\prime }=q(x,\theta ),(x,\theta )\in \Omega \times {\mathcal{S}}^{n-1},$$

$$f(\theta ,{\theta }^{\prime })=\{\begin{array}{rlrl}& \frac{1-g^2}{2\pi (1+g^2-2g\theta \cdot {\theta }^{\prime })},& & n=2,\\ & \frac{1-g^2}{4\pi (1+g^2-2g\theta \cdot {\theta }^{\prime }{)}^{3/2}},& & n=3,\end{array}$$

$$(\mathbf{K}\varphi )(x,\theta ):={\oint }_{{\mathcal{S}}^{n-1}}f(\theta ,{\theta }^{\prime })\varphi (x,{\theta }^{\prime })\mathrm{d}{\theta }^{\prime }.$$

$$\varphi (x,\theta )=q_b(x,\theta ),(x,\theta )\in {\Gamma }_{-},$$

$$h(x):={\mu }_a(x)\Phi (x),$$

$$\Phi (x):=(\mathbf{A}\varphi )(x)={\oint }_{{\mathcal{S}}^{n-1}}\varphi (x,\theta )\mathrm{d}\theta ,$$

$$p_0(x):=\gamma (x)h(x),$$

$$\{\begin{array}{rlrl}& \frac{1}{c^2(x)}\frac{{\partial }^2}{\partial t^2}p(x,t)-\Delta p(x,t)=0,& & (x,t)\in {\mathcal{R}}^n\times (0,T],\\ & p(x,0)=p_0(x),& & x\in \Omega ,\\ & \frac{\partial p}{\partial t}(x,0)=0,& & x\in \Omega ,\end{array}$$

$$\{\begin{array}{rlrl}& \left[\theta \cdot \nabla +{\mu }_a(x)+{\mu }_s(x)\right]\varphi (x,\theta )-{\mu }_s(x)(\mathbf{K}\varphi )(x,\theta )=0,& & (x,\theta )\in \Omega \times {\mathcal{S}}^{n-1},\\ & \varphi (x,\theta )=q_b(x,\theta ),& & (x,\theta )\in {\Gamma }_{-},\end{array}$$

$$\{\begin{array}{rlrl}& \left[\theta \cdot \nabla +{\mu }_a^2+{\mu }_s\right]\varphi (x,\theta )-{\mu }_s(\mathbf{K}\varphi )(x,\theta )=({\mu }_a^2-{\mu }_a^1){\varphi }^1(x,\theta ),& & (x,\theta )\in \Omega \times {\mathcal{S}}^{n-1}\\ & \varphi (x,\theta )=0,& & (x,\theta )\in {\Gamma }_{-}\end{array}$$

$${\varphi }^1\ge {\varphi }^2$$

$${\mu }_a^0(x)\le {\mu }_a^1(x)\le \cdots \le {\mu }_a^i(x)\le \cdots \le {\mu }_a^{\ast }(x).$$

$${\mu }_a^0(x)\le {\mu }_a^1(x)\le \cdots \le {\mu }_a^i(x)\le {\mu }_a^{\ast }(x).$$

$${\widetilde{{\mu }_a}}^{i+1}(x)\le {\mu }_a^{\ast }(x).$$

$${\mu }_a^i(x)\le {\mu }_a^{i+1}(x)\le {\mu }_a^{\ast }(x).$$

$${\varphi }^0(x,\theta )\ge {\varphi }^1(x,\theta )\ge \cdots \ge {\varphi }^i(x,\theta )\ge \cdots \ge {\varphi }^{\ast }(x,\theta ).$$

$$\{\begin{array}{rl}{\Omega }_{+}^i& :=\{x\in \Omega :h^{\ast }(x)\ge h^i(x)\},\\ {\Omega }_{-}^i& :=\{x\in \Omega :h^{\ast }(x)<h^i(x)\}.\end{array}$$

$${\widetilde{{\mu }_a}}^{i+1}$$

$$=\{\begin{array}{rlrl}& \ge {\mu }_a^i,& & x\in {\Omega }_{+}^i,\\ & <{\mu }_a^i,& & x\in {\Omega }_{-}^i.\end{array}$$

$${\mu }_a^{i+1}(x)=\{\begin{array}{rlrl}& {\widetilde{{\mu }_a}}^{i+1}(x),& & x\in {\Omega }_{+}^i,\\ & {\mu }_a^i,& & x\in {\Omega }_{-}^i.\end{array}$$

$$\underset{i\to \infty }{\lim }\Vert h^{\ast }-h^i{\Vert }_1=0.$$

$$h^{I+1}(x)={\mu }_a^{I+1}(\mathbf{A}{\varphi }^{I+1})(x)\le {\mu }_a^{I+1}(x)(\mathbf{A}{\varphi }^I)(x)=h^{\ast }(x).$$

$${\mu }_a^{i+1}(x)-{\mu }_a^i(x)$$

$$=\frac{h^{\ast }(x)-h^i(x)}{\mathbf{A}{\varphi }^i(x)}$$

$$\ge \frac{h^{\ast }(x)-h^i(x)}{2\pi \overline{{\varphi }^0}},$$

$${\varphi }^i(x)\le {\varphi }^0(x)\le \overline{{\varphi }^0}$$

$$\overline{{\mu }_a}(x)-{\mu }_a^0(x)=\sum _{i=0}^{\infty }({\mu }_a^{i+1}(x)-{\mu }_a^i(x))\ge \frac{1}{2\pi \overline{{\varphi }^0}}\sum _{i=0}^{\infty }(h^{\ast }(x)-h^i(x)).$$

$$h^{\ast }(x)-h^i(x)\stackrel{i\to \infty }{\to }0(\forall x\in {\Omega }_{+}^0),$$

$$\Vert {\mu }_a^{i+1}-{\mu }_a^i{\Vert }_{L^1({\Omega }_{+}^0)}\ge \frac{1}{2\pi \overline{{\varphi }^0}}\Vert h^{\ast }-h^i{\Vert }_{L^1({\Omega }_{+}^0)}.$$

$$\Vert \overline{{\mu }_a}-{\mu }_a^0{\Vert }_{L^1({\Omega }_{+}^0)}$$

Full text

1 Introduction

Photoacoustic tomography (PAT) is a developing medical imaging technique using non-ionizing waves. It combines ultrasound imaging and optical imaging [29], so it can achieve high contrast and high resolution simultaneously. Photoacoustic imaging is based on photoacoustic effect. A short pulse of electromagnetic energy propagates through the tissue rapidly at speed of higher magnitude than sound after the tissue is illuminated. Then the electromagnetic energy is absorbed by the tissue and converted to heat. There quickly follows the temperature of tissue increasing which leads to expansion spatially, and the expansion induces a pressure field. Acoustic waves generated by pressure distribution will be measured by detectors placed at the surface. They can be used to extract physical information about tissue inside. In medical imaging science, this technique is mostly applied to tomography imaging of skin, breast cancer and small animal imaging [29].

In the physical view, there are two inverse problems associated with two stages: the first stage of generating pressure after absorbing electromagnetic energy and second stage of initial pressure generating acoustic waves. First problem (aka PAT) is to reconstruct initial pressure distribution from temporal measurements of photoacoustic waves, and the second one, called quantitative photoacoustic tomography (QPAT), is to estimate optical coefficients from initial pressure distribution. Although initial pressure reflects some information inside, it is the product of the absorption of energy relied on properties of the object indirectly. Only the intrinsic optical coefficients provide direct information. Essentially, spatially varying optical coefficients lead to spatially varying initial pressure distribution. Moreover, reference [22] mentions that initial pressure decreases quickly in the tissue area with small absorption coefficient, so that we can not distinguish visually in the finial image of initial pressure. In a word, the original information in the investigated object can not be obtained completely until the optical coefficients are recovered. Generally, estimated optical coefficients comprise of absorption coefficient, scattering coefficient and Grüeneisen coefficient, where absorption coefficient is usually of the major clinic interest [14]. Nowadays, algorithms on reconstruction of initial pressure have been developed maturely [20, 29, 28, 1, 21], so a lot of researchers are focusing on QPAT.

In terms of mathematical model of QPAT, light propagation inside the tissue can be modeled by radiative transfer equation (RTE) [25, 22, 18] or diffusion approximation (DA) [4, 12]. However, DA approximates to RTE only in the situation where light behaves diffusely, which is invalid for small object or in regions close to light source, so that the optical coefficients estimated from DA model are not accurate. The reconstruction theories and some algorithms based on DA model are presented in [23, 24, 5, 13, 31, 7]. The two models are compared numerically in [27] and the derivation from RTE model to DA model can be found in [10].

In terms of measurement, there are mainly three kinds of methodologies: single illumination, multi-source illumination and multi-spectral illumination. It has been proved that single illumination can recover any one of absorption and Grüeneisen coefficient stably given scattering coefficient, and even there is an explicit analytic formula for the non-scattering media [22]. Multi-source can recover two of the three coefficients mentioned above uniquely and simultaneously [22, 18, 3, 26]. Furthermore, multi-spectral illumination can achieve unique reconstruction of three coefficients [2]. We focus on the RTE model with multi-source illumination and aim to investigate more efficient reconstruction method in this setting.

In term of reconstruction algorithm, assuming the knowledge of scattering coefficient, fixed-point iteration is used to recover absorption coefficient in [19]. However, the requirement of close enough initial value and the lack of theoretical results are obstacles of its application. When the scattering coefficient is unknown, QPAT can be formulated as a nonlinear optimization problem [18, 14] and furthermore can be simplified to a linear optimization problem by Born linearization [22]. Optimization approaches such as Jacobian-based method and gradient-based method are investigated in reference [25], in which Limited-memory BFGS (LBFGS), as a state of art gradient-based method, has been applied to the problem [14].

In the optimization framework, it is inevitable to solve original RTE in forward problem and adjoint RTE in the process of calculating gradient of objective function. Finite element method combined with streamline diffusion modification is applied to solve stationary RTE [18, 25, 30]. In this way, a large sparse linear system needs to be solved, which is still difficult. Inspired by [15], we apply Discontinuous Galerkin (DG) method combined with multigrid method to solve the two RTEs.

Our contributions in this paper are in several folds.

We propose an improved fixed-point iterative algorithm and prove its convergence given scattering coefficient. Numerically, the improved fixed-point iteration converges faster than the optimization-based Barzilai-Borwein (BB) method. 2. 2.

We apply BB method to reconstruct absorption coefficient and scattering coefficient simultaneously. To the best of our knowledge, it is the first time to apply BB method to this problem. 3. 3.

We revise the numerical scheme for RTE proposed in [15], which fails to solve adjoint RTE. We revise it to make it suitable for the adjoint RTE and prove its convergence. Then the RTE can be reduced to a sparse block diagonal linear system finally.

The rest of the paper is organized as follows. We introduce mathematical model of forward problem and inverse problem as well as notation conventions in Section 2. Given scattering coefficient, an improved fixed-point iterative algorithm and relevant proofs for convergence are presented in Section 3.1. The BB method for QPAT and the deduction of the gradient of objective function are detailed in Section 3.2. And we discuss solvers for original RTE and adjoint RTE respectively in Section 4.1. The simulation results are shown in Section 4.2. The conclusion comes in Section 5.

2 The mathematical description

In the following, we consider the convex bounded object region $\Omega\in\mathcal{R}^{n}(n=2,3)$ with Lipschitz boundary $\partial\Omega$ and angular space $\mathcal{S}^{n-1}$. Describing boundary condition, we divide boundary $\partial\Omega$ into inflow boundary and outflow boundary

[TABLE]

where $\nu(x)$ is the outward normal vector at the position $x$ of the boundary $\partial\Omega$.

2.1 QPAT model

We use stationary RTE to describe the light propagation in the form of

[TABLE]

where $\theta$ and $x$ are the direction and the position of interest respectively, $\phi(x,\theta)$ is the density of energy at the position $x$ in the direction $\theta$, $\mu_{a}(x)$ and $\mu_{s}(x)$ are spatially varying absorption and scattering coefficients respectively, and $q(x,\theta)$ is the source term. Function $f(\theta,\theta^{\prime})$ describing the probability of photon traveling from direction $\theta^{\prime}$ to $\theta$ is called scattering phase function and usually is characterized by anisotropic factor $g$ of the form

[TABLE]

which is the well-known Henyey-Greenstein (H-G) scattering function. For the sake of simplicity, we define scattering operator $\bm{K}$ by

[TABLE]

We assume inflow boundary condition

[TABLE]

that is no photon getting in $\Omega$ except for boundary source $q_{b}(x,\theta)$. In QPAT, source term $q(x,\theta)$ is usually regarded as zero. Absorbed energy density is presented by

[TABLE]

where

[TABLE]

in which $\bm{A}$ is the accumulation operator over all directions.

Owing to absorbed energy, initial acoustic pressure distribution is generated of the form

[TABLE]

where $\gamma(x)$ is spatially varying Grüeneisen coefficient. $\gamma(x)$ is the efficiency of conversion of absorbed energy to pressure and is assumed to be known. Throughout this work, it is rescaled to be one, i.e. $\gamma(x)=1$.

2.2 PAT model

On account of initial pressure distribution, acoustic waves propagate through object region $\Omega$ and a series of temporal acoustic signals are detected by ultrasonic detectors at the surface $\partial\Omega$. The behavior is described by

[TABLE]

where $c(x)$ is the acoustic speed, $p(x,t)$ is the acoustic pressure at the position $x$ and time $t$, and $p_{0}(x)$ is the initial pressure distribution.

The goal of PAT is to recover optical coefficients $\mu_{a}$ and $\mu_{s}$ given temporal data $p(x,t)\ ((x,t)\in\partial\Omega\times(0,T])$. This is often accomplished in two stages: first one recovers $p_{0}(x)$ from temporal data $p(x,t)\ ((x,t)\in\partial\Omega\times(0,T])$, and next one reconstructs $\mu_{a}$ and $\mu_{s}$ from initial pressure $p_{0}$ (or $h$) given illumination conditions. We assume the first stage has been done, and the second stage QPAT is our concern. We assume tissue is illuminated for $M$ times with known boundary sources $q_{b_{m}}\ (m=0,1,\dots,M-1)$, then inverse problem QPAT is to determine $\mu_{a}$ and $\mu_{s}$ from data $h_{m}$, that is $h(x;\mu_{a},\mu_{s},q_{b_{m}})\ (m=0,1,\dots,M-1)$.

2.3 Notation conventions

We denote the absorption coefficient and scattering coefficient by $\mu_{a}^{i}$ and $\mu_{s}^{i}$ in the $i$th iteration respectively, and corresponding solution of RTE and heat with $m$th illumination of boundary source $q_{b_{m}}\ (m=0,1,\dots,M-1)$ is denoted by $\phi_{m}^{i}(x,\theta)$ and $h^{i}_{m}(x)$. The exact coefficients and data are denoted by $\mu_{a}^{*}$, $\mu_{s}^{*}$, $\phi^{*}_{m}$, and $h^{*}_{m}$.

3 Reconstruction for optical coefficients

In this section, we propose an improved fixed-point iterative method to recover absorption coefficient given scattering coefficient and then apply Barzilai-Borwein (BB) method to simultaneously recover scattering coefficient as well. First, it is necessary to make several assumptions on the two optical coefficients.

Assumption 1**.**

Assume optical coefficients $\mu_{a}$, $\mu_{s}$ and energy density $\phi_{m}(x,\theta;\mu_{a},\mu_{s})$, which is the solution of RTE with $\mu_{a}$ and $\mu_{s}$, satisfy

$\mu_{a}(x)\in\mathcal{D}_{a}:=\{\mu_{a}\in C(\overline{\Omega}):0<\mu_{a}^{0}\leq\mu_{a}(x)\leq\mu_{a}^{\text{upper}}\}$* for fixed $\mu_{a}^{\text{upper}}$;* 2. 2.

$\mu_{s}(x)\in\mathcal{D}_{s}:=\{\mu_{s}\in C(\overline{\Omega}):0<\mu_{s}^{0}\leq\mu_{s}(x)\leq\mu_{s}^{\text{upper}}\}$* for fixed $\mu_{s}^{\text{upper}}$;* 3. 3.

Boundary condition $q_{b_{m}}(m=0,1,\dots,M-1)\in L^{\infty}(\Gamma_{-},|\theta\cdot\nu|)$; 4. 4.

Scattering kernel $f(\theta,\theta^{\prime})\in L^{1}(\mathcal{S}^{n-1}\times\mathcal{S}^{n-1})$; 5. 5.

If $\mu_{a}$ and $\mu_{s}$ satisfy (i) and (ii), there exist boundary sources $q_{b_{m}}(m=0,1,\dots,M-1)$ such that $\phi_{m}(x,\theta)\geq\epsilon>0$ for some positive $\epsilon$.

3.1 The case of given scattering coefficient

Under the Assumption 1, we provide an effective fixed-point iterative method. Since one measurement can recover absorption coefficient, subscripts ‘$m$’ of symbols mentioned in Section 2.3 are omitted for brevity. The algorithm is detailed as follows.

Remark*.*

The Assumption 1 is important because $\phi$ is a denominator in the Algorithm 1. In practice, we can also use $\phi+\varepsilon$ instead with the small $\varepsilon$ to avoid the situation where $\phi=0$.

Fixed-point iterative method has been mentioned in many references [19, 9, 14, 11], which is usually unstable and sensitive to initial value. However, the modification of fifth step in Algorithm 1 can guarantee convergence given scattering coefficient.

Carrying out Algorithm 1, we can obtain function sequences $\{\widetilde{\mu_{a}}^{i}(x)\}_{i=0}^{\infty}$, $\{\mu_{a}^{i}(x)\}_{i=0}^{\infty}$, $\{h^{i}(x)\}_{i=0}^{\infty}$ and $\{\phi^{i}(x)\}_{i=0}^{\infty}$. We claim that the sequence $\left\lVert h^{i}-h^{*}\right\rVert_{1}$ converges to zero. The proof is presented as follows.

Firstly, we state a theorem on stationary RTE.

Theorem 1** ([8]).**

For RTE (2) satisfying Assumption 1, if source $q\geq 0$ and boundary source $q_{b}\geq 0$, there exists a unique non-negative and continuous solution $\phi(x,\theta)$.

In the view of physics, the intuitive explanation behind Theorem 1 is straight-forward. For non-negative initial energy, the energy will not become negative while it was transferred and absorbed. Lemma 2 reveals the monotonicity of $\phi$ with respect to $\mu_{a}$.

Lemma 2**.**

Let $\phi^{1}(x)$ and $\phi^{2}(x)$ be the solutions of RTEs

[TABLE]

with $\mu_{a}$ being $\mu_{a}^{1}$ and $\mu_{a}^{2}$ respectively. Then $\phi^{1}(x,\theta)\geq\phi^{2}(x,\theta)$ provided $\mu_{a}^{1}(x)\leq\mu_{a}^{2}(x)(\forall x\in\Omega)$. Note that the superscripts of $\mu_{a}^{1}$ and $\mu_{a}^{2}$ are used to distinguish the different absorption coefficients.

Proof.

We see that $\phi:=\phi^{1}-\phi^{2}$ solves the RTE

[TABLE]

with absorption coefficient $\mu_{a}^{2}$ and scattering coefficient $\mu_{s}$. Since source term $(\mu_{a}^{2}-\mu_{a}^{1})\phi^{1}$ and boundary source are both non-negative, we can derive

[TABLE]

using Theorem 1. ∎

Using Lemma 2, we investigate the properties of function sequences from Algorithm 1.

Lemma 3**.**

For any $x\in\Omega$, the sequence $\{\mu_{a}^{i}(x)\}_{i=0}^{\infty}$ obtained in the Algorithm 1 satisfies

[TABLE]

Proof.

Given $\mu_{a}^{0}(x):=\mu_{a}^{0}\leq\mu_{a}^{*}(x)$, let us assume

[TABLE]

Then it suffices to establish (13) for $(i+1)$.

Obviously, $\mu_{a}^{i}(x)\leq\mu_{a}^{*}(x)$ indicates $\phi^{i}(x)\geq\phi^{*}(x)$. Since $\widetilde{\mu_{a}}^{i+1}(x)=\frac{h^{*}(x)}{(\bm{A}\phi^{i})(x)}$ and $\phi^{i}\geq\phi^{*}$, combining $h^{*}=\mu_{a}^{*}(\bm{A}\phi^{*})$, we infer

[TABLE]

According to Algorithm 1,

[TABLE]

An easy induction completes the proof. ∎

Monotonicity of $\mu_{a}^{i}$ can easily deduce the monotonicity of $\phi^{i}$.

Lemma 4**.**

For any $(x,\theta)\in\Omega\times\mathcal{S}^{n-1}$, the sequence $\{\phi^{i}(x,\theta)\}_{i=0}^{\infty}$ obtained in the Algorithm 1 satisfies

[TABLE]

Proof.

It is obvious using Lemma 2 and Lemma 3. ∎

We analyze the update process for $\mu_{a}$ by introducing new definitions. For any integer $i$, we divide region $\Omega$ into two parts in the form of

[TABLE]

Then we have

[TABLE]

According Algorithm 1,

[TABLE]

Actually, it is a process of keeping and updating from $\mu_{a}^{i}$ to $\mu_{a}^{i+1}$, keeping the value of $\mu_{a}^{i}$ of $\Omega^{i}_{-}$, updating the value of $\Omega^{i}_{+}$.

Eventually, we prove the convergence of function sequence $\{h^{i}(x)\}_{i=0}^{\infty}$.

Theorem 5**.**

For any $x\in\Omega$, the sequence $\{h^{i}(x)\}_{i=0}^{\infty}$ obtained in the Algorithm 1 satisfies

[TABLE]

Proof.

We divide our proof into four claims.

(i): Function sequences $\{\mu_{a}^{i}(x)\}_{i=0}^{\infty}$, $\{\phi^{i}(x,\theta)\}_{i=0}^{\infty}$ and $\{h^{i}(x)\}_{i=0}^{\infty}$ converge pointwisely.

Since monotonicity and boundedness of $\{\mu_{a}^{i}(x)\}_{i=0}^{\infty}$ and $\{\phi^{i}(x,\theta)\}_{i=0}^{\infty}$ for any $(x,\theta)\in\Omega\times\mathcal{S}^{n-1}$, it follows that they converge pointwisely, whose limits are denoted by $\overline{\mu_{a}}(x)$ and $\overline{\phi}(x,\theta)$ respectively. Obviously, $\{h^{i}(x)\}_{i=0}^{\infty}$ converges pointwisely and its limit is denoted by $\overline{h}(x)$. And obviously, $\overline{h}(x)=\overline{\mu_{a}}(x)(\bm{A}\overline{\phi})(x)$.

(ii): For any $x\in\Omega$, if there is some $I$ such that if $x\in\Omega^{I}_{+}$, then $x\in\Omega^{i}_{+}(i>I)$.

It is obvious that $\mu_{a}^{I+1}(x)=\widetilde{\mu_{a}}^{I+1}(x)\geq\mu_{a}^{I}(x)$. According to Lemma 4, we can obtain

[TABLE]

And an easy induction shows that $x\in\Omega^{i}_{+}(i>I)$.

(iii): $\overline{h}(x)=h^{*}(x)\ (x\in\Omega_{+}^{0})$, $\|h^{*}-h^{i}\|_{L^{1}(\Omega_{+}^{0})}\rightarrow 0\ (i\rightarrow\infty)$.

Thanks to $x\in\Omega^{i}_{+}(i>0)$ for any $x\in\Omega^{0}_{+}$, it follows that

[TABLE]

where $\overline{\phi^{0}}=\sup_{x\in\Omega}\phi^{0}(x)$. According to Lemma 4,

[TABLE]

gives last inequality in (16). Then

[TABLE]

Owing to $0\leq\overline{\mu_{a}}-\mu_{a}^{0}<\infty$, we can obtain that

[TABLE]

that is, $h^{*}(x)-h^{i}(x)$ converges to 0 pointwisely. So $\overline{h}(x)=h^{*}(x)\ (x\in\Omega_{+}^{0})$. Integrating (16) with respect $x$ over $\Omega_{+}^{0}$ gives

[TABLE]

Since

[TABLE]

and $0\leq\|\overline{\mu_{a}}-\mu_{a}^{0}\|_{L^{1}(\Omega_{+}^{0})}<\infty$, directly,

[TABLE]

(iv): $\|h^{*}-h^{i}\|_{L^{1}(\Omega_{-}^{0})}\rightarrow 0\ (i\rightarrow\infty$).

We divide region $\Omega_{-}^{0}$ into two parts of form

[TABLE]

According the proof of (iii) and the definition of $\Omega^{i}_{-}$, it can easily be seen that

[TABLE]

and

[TABLE]

From (19), it is clear that

[TABLE]

Considering RTEs with absorption coefficient $\mu_{a}^{i}$ and $\mu_{a}^{*}$, it is easy to obtain

[TABLE]

Since phase function $f(\theta,\theta^{\prime})$ satisfies

[TABLE]

it follows that

[TABLE]

Integrating (22) with respect to $\theta$ gives

[TABLE]

Integrating (23) with respect to $x$ over $\Omega$ gives

[TABLE]

Applying Green’s formula to the second term of (24), we have

[TABLE]

Apparently, we have

[TABLE]

Naturally, when $i$ is sufficiently large, we have

[TABLE]

Hence, when $i$ tends to infinity, we have

[TABLE]

Thus

[TABLE]

This completes the proof of (iv).

Therefore,

[TABLE]

∎

Given scattering coefficient, by Algorithm 1, we can get monotonically increasing sequence $\mu_{a}^{i}(x)$ and monotonically decreasing sequence $\phi^{i}(x,\theta)$, which converge pointwisely to $\bar{\mu_{a}}\leq\mu^{*}$ and $\bar{\phi}\geq\phi^{*}$ respectively. Sequence $h^{i}(x)$ converges to exact function $h^{*}(x)$ in the $L^{1}$-norm.

3.2 Reconstruction of $\mu_{a}$ and $\mu_{s}$ simultaneously

In practice, the scattering coefficient is also unknown, so the improved fixed-point iteration method is not applicable for the two unknown coefficients case. It is necessary to establish a more general method to recover two coefficients simultaneously. From now on, we follow optimization approach to estimate $\mu_{a}$ and $\mu_{s}$. First, we define error functional

[TABLE]

to measure the distance between measurement data $h_{m}^{*}$ and estimated data $h_{m}$, where $h_{m}(x;\mu_{a},\mu_{s})$ equals $\mu_{a}(\bm{A}\phi_{m})(x;\mu_{a},\mu_{s})$ for estimated $\mu_{a}$ and $\mu_{s}$. Then the reconstruction of QPAT can be reformulated as

[TABLE]

Remark*.*

We replace $\left\lVert h_{m}(\ \cdot\ ;\mu_{a},\mu_{s})-h_{m}^{*}\right\rVert_{2}^{2}/2$ by $\left\lVert\log(h_{m}(\ \cdot\ ;\mu_{a},\mu_{s}))-\log(h_{m}^{*})\right\rVert_{2}^{2}/2$ as data-fidelity term is due to that the latter boosts the convergence of minimization method according to [27]. We discuss its advantage and properties in Appendix C.

Remark*.*

According to inverse problem theory, regularization is useful for handling ill-conditioned problem. Usually, in image science, ill-conditionedness results in the edges of image blur and noise amplifying. Even though QPAT is a typical nonlinear problem, as pointed out in [2, 22], it is relatively stable with the help of multiple measurements. We focus on the minimization of data-fidelity term for the reason that multiple measurements alleviate the ill-posedness numerically. The algorithms on the minimization of objective function with regularization term, such as some a priori information, can be easily derived by our framework. We do not explore this in this paper.

Using the definition of Fréchet derivative $\nabla\mathcal{F}$ of $\mathcal{F}$ in feasible direction $(h_{\mu_{a}},h_{\mu_{s}})\in\mathcal{D}_{a}\times\mathcal{D}_{s}$, we can write

[TABLE]

The solution of RTE with coefficients $\mu_{a}$ and $\mu_{s}$ can be regarded as a function with respect to $\mu_{a}$ and $\mu_{s}$, that is, $\phi=\phi(x,\theta;\mu_{a},\mu_{s})$. Fortunately, the directional derivative of $\phi(x,\theta;\mu_{a},\mu_{s})$ with respect to $\mu_{a}$ and $\mu_{s}$ in any feasible direction exists according to [18]. Furthermore, the directional derivative of data-fidelity with respect to $\mu_{a}$ and $\mu_{s}$ can be expressed analytically by $\mu_{a}$, $\phi(x,\theta;\mu_{a},\mu_{s})$ and its solution of adjoint RTE. This concludes in Proposition 6.

Similar to the proof of Proposition 3.3 in [18], we can express the gradient of objective functional as follows.

Proposition 6**.**

For any pairs $(\mu_{a},\mu_{s})\in\mathcal{D}_{a}\times\mathcal{D}_{s}$ and feasible direction $(h_{\mu_{a}},h_{\mu_{s}})\in\mathcal{D}_{a}\times\mathcal{D}_{s}$, we have

[TABLE]

where $\phi_{m}^{*}$ solves following adjoint RTE

[TABLE]

Notice that $\bm{A}^{*}$ is the adjoint operator of $\bm{A}$ and

[TABLE]

Proof.

Following the proof in [18] (Proposition 3.3), we have

[TABLE]

Notice that $\phi_{m}^{\prime}(\ \cdot\ ;\mu_{a},\mu_{s})(h_{\mu_{a}},h_{\mu_{s}})$ is the directional derivative of $\phi_{m}(\ \cdot\ ;\mu_{a},\mu_{s})$ with respect to $(\mu_{a},\mu_{s})$ in direction $(h_{\mu_{a}},h_{\mu_{s}})$.

It completes the proof. ∎

After deducing the gradient of error functional (26), all we need to do is to find appropriate stepsize in the negative gradient direction to decrease the functional. Many optimization methods can achieve this goal, including steepest descent, Quasi-Newton, and so on. However, these gradient-based methods usually involve the linesearch process. A linesearch step needs to solve original RTE or adjoint RTE for several times. Since solving RTE dominates the computational cost, linesearch is computationally intensive. To mitigate the heavy computational cost, the well-known BB gradient method is applied to compute stepsize. It is derived from a two-point approximation to the scant equation underlying Quasi-Newton methods [6, 17]. And it is R-superlinearly convergent in the two-dimensional quadratic case [6]. Without loss of generality, we denote $\mu_{a}$ or $\mu_{s}$ by $\mu$ , then the update formula at $k$th step is

[TABLE]

Generally, there are two choices about the stepsize $\alpha_{k}$:

[TABLE]

and

[TABLE]

where $s_{k}=\mu^{k}-\mu^{k-1}$ and $y_{k}=\nabla\mathcal{F}_{k}-\nabla\mathcal{F}_{k-1}$, and $\mathcal{F}_{k}$ is the iterative sequence of error functional $\mathcal{F}$, see (26). BB algorithm is detailed in Algorithm 2.

4 Numerical simulations

The reconstructions of absorption and scattering coefficients are investigated with simulations in two cases: the reconstruction of $\mu_{a}$ given $\mu_{s}$ and the reconstruction of $\mu_{a}$ and $\mu_{s}$ simultaneously.

4.1 Solver for RTE

We consider the numerical solver for RTE in 2-D. And it can be extended to 3-D with little effort. Finite element method combined with streamline diffusion modification is applied to solve stationary RTE [30, 18, 25], where $P_{1}$ Lagrangian elements in spatial and angular space are used. In this way, a large sparse linear system needs to be solved which is still difficult. By improving the algorithm proposed in [15], we can solve original RTE as well as adjoint RTE on 2D and 3D unstructured mesh by Discontinuous Galerkin (DG) method combined with multigrid method, which reduces the problem to a sparse block diagonal linear system.

We divide angular space into $P$ equal intervals, and the corresponding directions and angles are denoted by $\theta_{0},\theta_{1},\dots,\theta_{P-1}$ and $\beta_{0},\beta_{1},\dots,\beta_{P-1}$ with interval $\Delta\beta$. The spatial domain is discretized into unstructured triangular mesh. Suppose spatial domain is divided into N triangles and each triangle contains $n_{d}$ nodes, where $n_{d}=3$ for 2-D spatial domain. Lagrangian elements and piecewise linear DG elements are used to discretize RTE, that is, numerical discrete scheme of RTE is

where $L_{k}(\theta)$ is piecewise linear basis function in angular space which takes value of 1 in direction $\theta_{k}$ and 0 in other directions, $\varphi_{ij}$ is the linear basis function in spatial domain in the direction $\theta_{k}$ which takes value of 1 in the $j$th nodes of $i$th triangle and 0 in other nodes and triangles, and $\phi_{i,j,k}$ is the value of $\phi(x,\theta)$ in direction $\theta_{k}$ in $j$th node of $i$th element. Such spatial basis function can approximate discontinuous solution which is more suitable for actual situation, such as some edges of inclusions in object region. On account of scattering term in RTE, angular Gauss-Seidel iteration is applied in [15]. It iteratively solves RTE in every fixed direction in the form of

[TABLE]

where $\omega_{kk^{\prime}}=\frac{\omega^{0}_{kk^{\prime}}}{\sum_{k^{\prime}}\omega^{0}_{kk^{\prime}}}$ with $\omega^{0}_{kk^{\prime}}=f(\theta_{k},\theta_{k^{\prime}})$. Obviously, equation (33) can be solved by a lot of numerical methods, and we use DG method.

A multigrid scheme is applied to solve (33), and it is reduced to solve a sparse block diagonal system. Multiplying test function $\varphi_{ij}$ in both sides of (33), integrating it over $i$th triangle $\tau_{i}$ with respect to $x$, we have

[TABLE]

from Green’s formula, where subscripts $i$ and $j$ of $\phi_{i,j,k}$ are omitted. In (34), $\hat{\phi}$ is the value of neighboring element in inflow direction, which is the product of upwind scheme. The specific explanation of why the value inflow direction is used is as follows.

For an intuitive explanation, we use rectangular mesh to illustrate in Figure 1. For rectangular mesh, in direction $\theta_{k}$ and node $x$, RTE can be discretized into

[TABLE]

where $\phi_{k}:=\phi(x,\theta_{k})$. Backward difference is used to approximate one order derivative. Then when $0\leq\beta_{k}<\frac{\pi}{2}$, its upwind scheme is

[TABLE]

where

[TABLE]

Obviously, the scheme converges because the equation (36) is diagonally dominant, and Gauss-Seidel scheme will accelerate the convergence. In fact, from (36), $\phi_{i,j,k}$ is updated using

[TABLE]

Apparently, $\phi_{0,j,k}$ and $\phi_{i,0,k}$ are known as inflow boundary condition. Therefore, (37) is an explicit scheme which updates $\phi_{i,j,k}$, $\phi_{i,1,k}$, $\phi_{2,j,k}$, $\phi_{i,2,k},\cdots\ (i,j,=1,2,\cdots)$ successively. We can see that (36) and (34) both use inflow information to update outflow information. This may be a kind of intuitive explanation. Indeed, the convergence of (34) for vacuum boundary condition is proved in [16], see Appendix A.

However, the above numerical scheme fails to solve the adjoint RTE with same updating order, that is to use inflow information to update outflow information. We provide a reverse updating order for elements to guarantee the convergence of solver for adjoint RTE. The solver for adjoint RTE is detailed as follows.

For adjoint RTE (30), the corresponding discrete scheme in direction $\theta_{k}$ for rectangular mesh is

[TABLE]

Differently, forward difference is used to approximate the one order derivative, and when $0\leq\beta_{k}<\frac{\pi}{2}$, its upwind scheme is

[TABLE]

where

[TABLE]

Then $\phi_{i,j,k}$ is updated by

[TABLE]

The use of forward difference makes (39) diagonally dominant again. Besides, outflow information is used to update inflow information in (40). Heuristically, for unstructured mesh, we have

[TABLE]

where $\hat{\phi}$ is the value of neighboring element in outflow direction which is used to update inflow information. Using similar discussion, the convergence of (41) is proved in Appendix B. Furthermore, we apply multigrid method to solve RTE to accelerate convergence.

From above discussion, no matter which mesh is used, the key of convergence is to update outflow information using inflow information for the original RTE (2) and update the latter using the former for the adjoint RTE (41). The convergence in rectangular mesh is obvious. As for triangular mesh, the update scheme results from variation analysis, so the convergence proofs are obtained by discussing corresponding variations, see Appendix A and B.

Therefore, for original RTE (2), there are three layers of loops. We apply the multigrid scheme at the outermost loop. Since it involves two coordinate systems, i.e., angular-coordinate and spatial-coordinate, there are variable updating schemes in terms of the iteration order, such as the angle-prior and space-prior. The second-layer loop is about direction, that is, the corresponding spatial equation (33) is solved in turn for each direction $\theta_{k}(k=0,1,\dots,P-1)$. In third-layer loop, for each element in the direction $\theta_{k}$, the value of $\phi_{i,j,k}$ is updated iteratively through solving a $3\times 3$ linear system (34). Note that the updating order is in the following order: consider the each element successively from the boundary along the direction of $\theta_{k}$. We refer the interested reader to [15] for its algorithm details and to [16] for its theory. On the contrast, for adjoint RTE, in third-layer loop, inspired by (39) we propose to consider each element successively from the boundary along the direction of $-\theta_{k}$, which just reverse the updating order of original RTE. And linear system (41) needs to be solved. The corresponding pseudo-code is presented in Algorithm 3.

4.2 Numerical results

In this subsection, we present some numerical results to demonstrate the numerical performance of improved fixed-point iteration and BB method. For the sake of simplicity, only 2D reconstruction is investigated. The anisotropic factor $g$ equals 0.9. Applying RTE solver described in 4.1, we can get discrete energy density $\bm{H}$ of the same length as mesh. In order to explore the stability to noise, noisy data are generated by

[TABLE]

where $\bm{N}$ is a random vector that follows the normal distribution with mean 0 and variance 1. We consider 4 object regions:

Rectangle $\Omega_{0}=[-20,20]\times[-20,20]$ and four inclusions $\Omega_{1}=\{(x,y)\in\Omega_{0}:(x+10)^{2}+(y-10)^{2}=6^{2}\}$, $\Omega_{2}=\{(x,y)\in\Omega_{2}:(x-10)^{2}+(y-10)^{2}=4^{2}\}$, $\Omega_{3}=[-17,-5]\times[-17,-5]$, and $\Omega_{4}=[5,17]\times[-17,-5]$; 2. 2.

Circle $\Omega_{0}$ with center $(0,0)$ and radius 20; 3. 3.

Circle $\Omega_{0}$ with center $(0,0)$ and radius 20, and four inclusions $\Omega_{1}=[-12,-8]\times[-12,12]$, $\Omega_{2}=[-8,-2]\times[-12,12]$, $\Omega_{3}=[-2,12]\times[6,12]$, and $\Omega_{4}=[-2,12]\times[-12,6]$; 4. 4.

Circle $\Omega_{0}$ with center $(0,0)$ and radius 20, and three inclusions $\Omega_{1}=\{(x,y)\in\Omega_{0}:\frac{(y-3)^{2}}{9^{2}}+\frac{(x-7)^{2}}{6.2^{2}}=1\}$, $\Omega_{2}=[-14,-4]\times[-10,8]$, and $\Omega_{3}=\{(x,y)\in\Omega_{0}:\frac{(\frac{\sqrt{2}x}{2}+\frac{\sqrt{2}y}{2}-8.4)}{8^{2}}+\frac{(-\frac{\sqrt{2}x}{2}+\frac{\sqrt{2}y}{2}+8)^{2}}{5^{2}}=1\}$.

Their absorption and scattering coefficients respectively are

$\mu_{a}$ is 0.02 in $\Omega_{1}$ and $\Omega_{4}$ with background 0.01; $\mu_{s}$ is 3 in $\Omega_{2}$ and $\Omega_{3}$ with background 1, 2. 2.

$\mu_{a}(x,y)=0.02+0.01\sin(\frac{\pi x}{8})$ and $\mu_{s}(x,y)=2+\sin(\frac{\pi y}{8})$, 3. 3.

$\mu_{a}$ is 0.03 in $\Omega_{1}$, 0.02 in $\Omega_{2}$, 0.04 in $\Omega_{3}$, and 0.015 in $\Omega_{4}$ with background 0.01; $\mu_{s}$ is 2.5 in $\Omega_{1}$, 1.5 in $\Omega_{2}$, 3 in $\Omega_{3}$, and 2 in $\Omega_{4}$ with background 1, 4. 4.

$\mu_{a}$ is 0.015 in $(\Omega_{1}\cup\Omega_{2})\backslash(\Omega_{1}\cap\Omega_{2})\cup\Omega_{3}$ and 0.03 in $\Omega_{1}\cap\Omega_{2}$ with background 0.01; $\mu_{s}$ is 3 in $\Omega_{1}\cup\Omega_{3}$ with background 1.

These regions are depicted in Figure 2, where the optical coefficients of 1st, 3rd, and 4st templates are piecewise constant, and the second one is smooth in object region $\Omega_{0}$. These templates contain discontinuous borders as well as continuous borders, whose corners are both sharp and smooth. To avoid inverse crime, we generate data by solving RTE in finer unstructured mesh than inverse problem. Original data are generated on 21376, 16352, 17376, and 16576 unstructured mesh respectively, and corresponding inverse problem are solved on 9600, 7392, 7392, 7392 unstructured mesh. Four point sources are placed in (-20,0), (0,20), (20,0) and (0,-20). Iterative relative errors are defined by

[TABLE]

4.2.1 Reconstruction of $\mu_{a}$ given $\mu_{s}$

Using one measurement with boundary source in (20,0), we apply improved fixed-point iterative method and BB method to retrieve absorption coefficients of templates given scattering coefficients. Our initial guesses are set to be the same as background. Figure 3 shows results after 50 iterations for noiseless data. Then $5\%$ Gaussian noise (i.e. $\epsilon=5\%$) is added to data according to (42), and the final reconstructed images are showed in Figure 4. Specific iterative relative errors are showed in Figure 5. From Figure 3 and Figure 4, it seems that both methods can retrieve absorption coefficient with almost the same accurate solutions in the case of noise-free and noisy measurements. However, from Figure 5, it is obvious that improved fixed-point iteration converges more rapidly than BB method. Improved fixed-point iteration achieve the critical point only after a few iterations. To attain the same accuracy, the number of iteration of improved fixed-point iteration is about the half of the one of BB method. Owing to solving the adjoint RTEs in BB method, improved fixed-point iteration can achieve the accuracy as the BB method with about quarter computational time of the latter.

4.2.2 Reconstruction of $\mu_{a}$ and $\mu_{s}$ simultaneously

We apply BB method to reconstruct absorption and scattering coefficients simultaneously from four measurements. The boundary point sources are placed in the top, bottom, left and right sides respectively. The reconstruction results from noiseless data are showed in Figure 6. Corresponding reconstruction results from data added by 5% Gaussian noise are showed in Figure 7. The relative $L^{2}$ errors of reconstructions of optical coefficients of four templates are tabulated in Table 1. From Figure 6 and Figure 7, we can see that BB method retrieves absorption coefficient accurately in piecewise constant and smooth cases, and even is stable to noisy data. In the first, third and fourth templates, some jagged borders of inclusions are reasonable, which are caused by the interpolation of data from fine mesh to coarse mesh. Nevertheless, the reconstruction of scattering coefficient is not satisfactory. In piecewise constant case, the edges between different pieces are blur, which can be seen from the first, third, and fourth template in Figure 6 and Figure 7. This is due to the insensitive scattering coefficient. After all, there is no explicit $\mu_{s}$ in the formula (6).

5 Conclusion

In this paper, we investigate the reconstruction of absorption and scattering coefficients in QPAT in two cases. Given scattering coefficient, we propose an improved fixed-point iterative method to reconstruct the absorption coefficient and prove its convergence. The advantage of this reconstruction algorithm is that its fast convergence and it does not require the initial guess close to the exact solution. Meanwhile, it does not need to solve adjoint RTE in the optimization approach. For the simultaneous reconstruction of the two coefficients, we apply a state of art BB method, which dose not need linesearch to compute stepsize. Indeed, linesearch involves solving RTE for several times, which is expensive. BB method only use the values of previous two steps to update current estimation. Numerical results show that improved fixed-point iteration can achieve almost the same accuracy with BB method. It takes about quarter the computational time of the BB method. Moreover, the algorithm is stable to noise. In the unknown scattering coefficient case, numerical results show that BB method can obtain quite accurate absorption estimate. However, the result of the estimation of scattering coefficient is not satisfactory. There is no explicit $\mu_{s}$ in (6), so it is insensitive to the objective function. There is obvious blur near the borders in piecewise constant templates.

In future, we will design a better error functional to improve the reconstruction of $\mu_{s}$. Although BB method can recover numerically absorption and scattering coefficients simultaneously, its convergence remains to be studied. Besides, for more practical application, we will do some realistic simulations, such as 3-D templates or real physical QPAT data.

6 Acknowledgments

We would like to thank the anonymous referees for their useful comments that help us improve the quality of the paper. We would like to thank Prof. Markus Haltmeier (Department of Mathematics, University of Innsbruck Technikestraße) and Dr. Hao Gao (Wallace H. Coulter Department of Biomedical Engineering, Georgia Institute of Technology) for their helpful advice on RTE solver. We also thank Mr. Ji Li (School of Mathematical Sciences, Peking University) for helpful comments on this manuscript draft. This work was supported by NSF grants of China (61421062, 11471024).

Appendix A Convergence of the solver for original RTE

We use the norm $\|\phi\|:=\sqrt{\sum_{k=0}^{P-1}\omega_{k}\oint_{\Omega}\phi_{k}^{2}}$ to estimate the accuracy of numerical solution, where $\omega_{k}$ is defined by $\omega_{k}=\oint_{\mathcal{S}^{n-1}}L_{k}(\theta)\,\mathrm{d}\theta$ with Lagrangian function $L_{k}(\theta)$ in $\mathcal{S}^{n-1}$. We denote the solution of original RTE (2) and angular discretized equation (33) by $\phi(x,\theta)$ and $[\phi]:=(\phi_{k}(x))_{k=0}^{P-1}$ respectively, and assume the angular and spatial mesh size are $h_{a}$ and $h$ respectively. Then there are some convergence results.

Theorem 7** ([16]).**

Assume $\phi(x,\theta)\in C^{2}(\Omega\times\mathcal{S}^{n-1})$, with vaccum boundary condition and sufficiently fine angular mesh

[TABLE]

The convergence of spatial discrete scheme (34) is discussed next. Assume triangulation $T_{h}$ with $h=\sup_{S\in T_{h}}diam(S)$, $V_{h}^{d}:=\{v:v|_{S\in T_{h}}\in P^{d}(S),v|_{\Omega\backslash S}=0\}$, where $P^{d}(S)$ is the space of $d$-degree polynomials. Let $(u,v)_{S}:=\int_{S}uv\,\mathrm{d}x$, $\Gamma_{h}:=(\cup_{S\in T_{h}}\partial S)\backslash\partial\Omega$ and $\partial\Omega^{k}_{-(+)}:=\{x\in\partial\Omega:\nu\cdot\theta_{k}\leq(>)0\}$. Then the convergence result is as follows.

Theorem 8** ([16]).**

With vaccum boundary condition, if $\phi_{k}^{h}\in V^{d}_{h}$ satisfies

[TABLE]

where

[TABLE]

with $\phi_{k}^{h^{\pm}}=\lim_{\epsilon\rightarrow 0^{\pm}}\phi_{k}^{h}(x+\epsilon\theta_{k})$, $\left<u,v\right>_{\Gamma_{h}^{k}}=\int_{\Gamma_{h}}uv|\theta_{k}\cdot\nu|$, and $\left<u,v\right>_{\partial\Omega^{k}_{-}}=\int_{\partial\Omega^{k}_{-}}uv|\theta_{k}\cdot\nu|$, then $\phi^{h}:=(\phi_{k}^{h})\in(V_{h}^{d})^{P}$ satisfies

[TABLE]

where the $|\phi|^{2}_{d+1}:=\sum_{k=0}^{P-1}\omega_{k}|\phi_{k}|^{2}_{H^{d+1}}$.

Combining (43) and (45), when $d=1$ we can obtain the estimate

[TABLE]

Obviously, the corresponding update scheme of (44) is (34) when $d=1$. Therefore the convergence proof of (34) is done.

Appendix B Convergence of the solver for adjoint RTE

As for adjoint RTE (30), we denote the solution of it and its angular discretized equation

[TABLE]

by $\tilde{\phi}(x,\theta)$ and $[\tilde{\phi}]:=(\tilde{\phi}_{k}(x))_{k=0}^{P-1}$ respectively. Through similar discussion, we obtain following results.

Theorem 9**.**

Assume $\tilde{\phi}(x,\theta)\in\mathcal{C}^{2}(\Omega\times\mathcal{S}^{n-1})$, with vaccum boundary condition and sufficiently fine angular mesh

[TABLE]

Theorem 10**.**

With vaccum boundary condition, if $\tilde{\phi}_{k}^{h}\in V^{d}_{h}$ satisfies

[TABLE]

where

[TABLE]

then $\tilde{\phi^{h}}:=(\tilde{\phi}_{k}^{h})\in(V^{d}_{h})^{P}$ satisfies

[TABLE]

where the $|\tilde{\phi}|^{2}_{d+1}:=\sum_{k=0}^{P-1}\omega_{k}|\tilde{\phi_{k}}|^{2}_{H^{d+1}}$.

Combining (49) and (50), when $d=1$ we can obtain the estimate

[TABLE]

Obviously, the corresponding update scheme of (49) is (41) when $d=1$. Therefore the convergence proof of (41) is done.

Appendix C Comments on the log-type error functional

According to [27], taking log accelerates the convergence of minimization method significantly. It is due to taking log changes the shape of contours of error functional. From Figure 2 in [27], the contours of error functional become less narrow after taking log. In numerical simulations, we find log-type error functional decreasing more fast. Indeed, we can explain it by a one-dimension example. Let $f(x)$ be a continuous function on bounded region $[a,b]$, which satisfies $0<f(x)<1$. Given $h^{*}=f(x^{*})\ (\forall x\in[a,b])$, we can apply least square method to recover $x^{*}$. Define error functional

[TABLE]

Then the gradients are

[TABLE]

Theorem 11**.**

For error functionals defined in (52), we claim that there exists $\delta>0$ such that if $x\in(x^{*}-\delta,x^{*}+\delta)$, $|\nabla\mathcal{J}_{1}(x)|\leq|\nabla\mathcal{J}_{2}(x)|$.

Proof.

Assume there exists some $\delta_{0}$ such that if $x\in(x^{*}-\delta_{0},x^{*}+\delta_{0})\backslash\{x^{*}\}$, $f(x)\neq h^{*}$. Otherwise, $|\nabla\mathcal{J}_{1}|=|\nabla\mathcal{J}_{2}|$ in $(x^{*}-\delta_{0},x^{*}+\delta_{0})$.

Since $f(x)$ is continuous, $f(x)\rightarrow f(x^{*})\ (x\rightarrow x^{*})$, then

[TABLE]

Therefore, for arbitrary $\epsilon>0$, there exists some $\delta_{1}>0$ such that if $x\in(x^{*}-\delta_{1},x^{*}+\delta_{1})\backslash\{x^{*}\}$,

[TABLE]

Dividing (54) by $f(x^{*})$, we can obtain

[TABLE]

Let $0<\epsilon\leq 1-f(x^{*})$, we can get

[TABLE]

Clearly, it follows that

[TABLE]

Let $\delta=\min\{\delta_{0},\delta_{1}\}$, then $|\nabla\mathcal{J}_{1}|\leq|\nabla\mathcal{J}_{2}|$ in $(x^{*}-\delta,x^{*}+\delta)$. ∎

Theorem 12**.**

For error functionals defined in (52), we claim that there exists $\delta>0$ such that functional

[TABLE]

is monotonically increasing with respect to $\alpha(x):=|x-x^{*}|$ in domain $[x^{*}-\delta,x^{*}+\delta]$.

Proof.

Assume there exists $\delta_{0}$ such that $f(x)\neq h^{*}$ in $[x^{*}-\delta_{0},x^{*}+\delta_{0}]$. Since $f$ is continuous, there exists some $\delta_{1}<\delta_{0}$ such that $|f(x)-h^{*}|$ is increasing with respect to $\alpha(x)$, that is, in interval $[x^{*}-\delta_{1},x^{*}+\delta_{1}]$, the closer between $x$ and $x^{*}$ then the smaller $|f(x)-h^{*}|$.

For arbitrary $x=x^{*}+\Delta x\in[x^{*}-\delta_{1},x^{*}+\delta_{1}]$, let $f(x)=h^{*}+\Delta f$. Then we have

[TABLE]

Clearly, $\log(1+\frac{\Delta f}{h^{*}})+\Delta f$ and $\log(1+\frac{\Delta f}{h^{*}})-\Delta f$ are monotonically increasing with respect to $\Delta f$. Combing $\log(1+\frac{\Delta f}{h^{*}})+\Delta f=\log(1+\frac{\Delta f}{h^{*}})-\Delta f=0$ when $\Delta f=0$, we obtain $\mathcal{D}(x)$ is non-negative and monotonically increasing with respect to $|\Delta f|$. Therefore, $\mathcal{D}$ is monotonically increasing with respect to $\alpha(x)$.

Let $\delta=\min\{\delta_{0},\delta_{1}\}$, then this completes the proof. ∎

From Theorem 11 and Theorem 12, we see if $\sup f(x)<1$, taking log not only makes the error functional steeper near the minimizer for some fixed direction, but also the farther between $x$ and $x^{*}$, the change from $\mathcal{J}_{1}(x)$ to $\mathcal{J}_{2}(x)$ is more significant. Then, for multivariate function, taking log makes the contours not so narrow. Thence, log-type error functional improves the convergence of minimization method. As for the general case of $\sup f(x)>1$, we can replace $f(x)$ by $f(x)/c$ for some constant $c>\sup f(x)$.