{'fraction_non_alphanumeric': 0.05171849427168576, 'fraction_numerical': 0.020990180032733223, 'mean_word_length': 3.951580226904376, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical reconstruction of radiative sources in an absorbing and non-diffusing scattering medium in two dimensions

We consider the two dimensional quantitative imaging problem of recovering a radiative source inside an absorbing and scattering medium from knowledge of the outgoing radiation measured at the boundary. The medium has an anisotropic scattering property that is neither negligible nor large enough for the diffusion approximation to hold. We present the numerical realization of the authors' recently proposed reconstruction method. For scattering kernels of finite Fourier content in the angular variable, the solution is exact. The feasibility of the proposed algorithms is demonstrated in several numerical experiments, including simulated scenarios for parameters meaningful in optical molecular imaging.

The medium is characterized by the absorption coefficient µ a (x), scattering coefficient µ s (x), and scattering kernel p(x, ξ · ξ ), all of which are assumed known, real valued, non-negative, essentially bounded functions. For any x ∈ D and ξ, ξ ∈ S 1 , the scattering kernel (a conditional probability) satisfies S 1 p(x, ξ · ξ )dσ ξ = 1, for all x ∈ D. The total attenuation is defined as µ t = µ a + µ s . Generated by an unknown source q, in the steady state case, and assuming no incoming radiation from outside the domain, the density of particles I(x, ξ) at x ∈ D traveling in the direction ξ satisfies the problem ξ · ∇ x I(x, ξ) + µ t (x)I(x, ξ) The unknown source q is assumed square integrable and compactly supported in D. Under the above assumptions on µ a , µ s , p, and q, the forward problem (1) is well posed ( [10]), with a unique solution I in the space For well-posedness results under various other assumptions, see [1,8,9,22], and the generic result in [35].
For a given medium, i.e., µ a , µ s and p are known, we consider the following inverse source problem: Determine q in D from the measurement I(x, ξ) = I measure (x, ξ), (x, ξ) ∈ Γ + of the directional outflow at the boundary. In the presence of scattering, this problem is equivalent to inverting a smoothing operator, and therefore it is illposed.
When µ a = µ s = 0, this is the classical X-ray tomography problem of Radon [31], where q is to be recovered from its integrals along lines, a problem which is well understood both on its theoretical and numerical facets, see, e.g., [16,18,26] and references therein. For µ a = 0 but µ s = 0, this is the problem of inversion of the Attenuated X-ray transform in two dimensions, solved by different methods in [2], and [28]; see [4,6,19,23,24,27] for later approaches and some numerical implementation.
The inverse source problem in a scattering media considered here, i.e., µ s p = 0, has also been considered under various limiting constraints (e.g., [5,17,20,21,34]) with a most general result showing that the source is uniquely and stably determined by the outflow in [35]. However, the numerical solutions for the inverse source problem based on the above mentioned results have yet to be realized.
In the special case of a weakly (anisotropic) scattering media, where 1−p 1, one may recover the source by devising algorithms for the iterative method proposed in [5]. However, based on a perturbation argument of the non-scattering case in [28], the method does not extend to strongly anisotropic scattering media considered here. Moreover, even in the case of a weakly scattering media, the requirement of solving one forward problem (a computationally extensive procedure) at each iteration renders the method inefficient from an imaging perspective.
In here we present a numerical reconstruction method based on the authors' recent theoretical results in [14], and propose one algorithm to recover the radiative sources. As demonstrated in the numerical experiments below, the algorithms can handle the quantitative imaging of sources in a non-small scattering medium that is far from diffusive approximation, with applications to Optical Molecular Imaging [29,30].
Key to the reconstruction method is the realization that any finite Fourier content in the angular variable of the scattering kernel splits the problem into a nonscattering one and a boundary value problem for a finite elliptic system. The role of the finite Fourier content has been independently recognized in [25]. However, in general, the scattering kernel does have an infinite Fourier content. Such is the case in the numerical experiments below, where we work with the ubiquitous (two dimensional version of) the Henyey-Greenstein kernel (2) p(x, ξ · ξ ) = p(ξ · ξ ) = 1 2π In (2) the parameter 0 ≤ g ≤ 1 models a degree of anisotropy with g = 0 being the isotropic case, and g = 1 being the ballistic case. In the proposed reconstruction we use an approximate scattering kernel obtained by truncating the Fourier modes in the angular variable. The error estimate in the data due to such truncation (see Section 3 below) allows to interpret the reconstruction as a minimum residual solution, where for a given a priori noise level, the degree of truncation dictates the number of significant Fourier modes. Moreover, we devise a locally optimal criterion for the choice of the order of this truncation, which is independent of the unknown source.
The theoretical method originated in Bukhgeim's theory of A-analytic functions developed in [7] to treat the non-attenuating case and in [2] for the absorbing but non-scattering case, and extends the ideas in [32,33] to the scattering case. One of the numerical results in the first experiment below has been announced in [14].

Preliminaries
In this section we establish notation, while presenting the basic ideas of the reconstruction in the non-scattering case. The presentation follows authors' ideas in [14,32,33], while the notation is that used by practitioners in Optical Tomography, e.g. [3]. We also depart from the original notation in [7] (used so far in [14,32,33]), and work with the positive Fourier modes. This will allow for the natural indexing of sequences. For the analytical framework, which specifies the regularity of the coefficients and proves the appropriate convergence of the ensuing series we refer to [14]. We identify points x = (x 1 , x 2 ) ∈ R 2 with their complex representative z = x 1 + ix 2 ∈ C. The method considers the transport model (1) in the Fourier domain of the angular variable. For ξ = ξ(θ) = (cos θ, sin θ) ∈ S 1 , let be the Fourier series representation of I in the angular variable. From the analysis of forward problem (e.g., [10]), for a square integrable source and essentially bounded coefficients, the solution I of (1) is at least square integrable in the angular variable and thus the series (3) is summable in L 2 sense. Moreover, since I is real valued, I −m = I m and the angular dependence is completely determined by the sequence of its nonnegative Fourier modes D z → I 0 (z), I 1 (z), I 2 (z), ... .
The system (5) was originally introduced (in a context more general than needed here) in [7], and shown that their solutions satisfy a Cauchy like integral formula, where the interior values in D are recovered from their boundary values. More precisely, for each m ≥ 0, and each z ∈ D fixed, In the absorbing and non-scattering case (µ s = 0 but µ a ≥ 0) the system (4) becomes The system (7) has been studied first in [2] and shown that an integrating operator can be found to reduce it to (5). We briefly describe the explicit construction of the integrating factor introduced in [12] and its convolution form in [32] to be used in our reconstruction algorithms below. For brevity assume that µ a and µ s are extended by zero outside the domain.
s − t dt be the Hilbert transform, where ξ ⊥ ∈ S 1 denotes the counterclockwise rotation of ξ by π/2. Let us define While there are many possible such integrating factors to reduce the system (7) to (5), the key feature of the construction in [12] is the fact that all the negative modes of θ → h[µ t ](x, ξ(θ)) (and thus of e ±h ) vanish, as shown in [6,12,26].
} be the corresponding sequences of the Fourier modes, Then, as shown in [32, Lemma 4.1], if (7), then { J m } is a solution of (5), and conversely, if (5), then { I m } is a solution of (7). Moreover, the Cauchy problem for (5) (and thus for (7)) has at most one solution.

Reconstruction in the presence of scattering
We consider now the scattering case, when the scattering kernel p is of polynomial type in the angular variable, i.e., where M is the degree of the polynomial. We stress that no smallness is assumed on the Fourier modes p k , for k = 0, ..., M .
In this case the transport equation (1a) reduces to the system The basic idea in the reconstruction starts from the observation that the system (12c) is of the type (7). Therefore, via the integrating formulas (10) for m ≥ M , the problem reduces to finding { J (M ) m } m≥M solution of the Cauchy problem for the elliptic system (5). Moreover, from the data on the boundary we can recover for each m ≥ M , and by using (10) we find the boundary data Next we use the Cauchy-like integral formula (6) to recover the interior values { J m (z) ; z ∈ D } for m ≥ M . Finally, we use (11) to recover the interior values { I m (z) ; z ∈ D }, Recursively and in the decreasing order starting with the index m = M − 1 to m = 0, we solve the elliptic problems (12b) as follows. By applying ∂ to (12b), we are lead to solving while on the boundary  , we can now use (12a) to reconstruct the unknown source q in D.
Remark. In general, the scattering kernel p is not of polynomial type. An immediate application of the well posedness of the forward problem in L 2 (D × S 1 ) gives an error estimate in the measured outflow due to the approximation in the scattering kernel. Proposition 1. Let µ t , µ s , p be such that the forward problem (1) has a unique solution and letp be such that for some 0 < < 1, and letĨ be the unique solution for (1) corresponding top. Then, where C > 0 is a constant depending only on µ t , µ s , p and the domain.
Proof. It is easy to see that the difference of the two corresponding solutions h = I −Ĩ satisfy By interpreting the last term as a source and using the classical estimates in the forward model; see, e.g., [35, Theorem 2.1], we obtain For the two dimensional Henyey-Greenstein kernel (2), with 0 ≤ g < 1, considered in the numerical simulations and its M -th order truncation (15) p one obtains a refined estimate in terms of the anisotropic parameter g. Namely, where the second last equality uses Parseval's identity. For a given level of noise, a sufficiently large choice of M yields that the difference (14) between the exact data and the hypothetical data (which would be obtained had the scattering been of polynomial type) falls under the noise level. Therefore our reconstruction produces a source for which the corresponding boundary data is indistinguishable from the exact data within the level of noise. This is the most we can hope to reconstruct. It is worth noting that, for the scattering kernels of polynomial type, the method above does produce the exact solution in a stable manner, as shown by the authors in [14, Corollary 6.1].

Evaluation of the Hilbert
Transform. An accurate calculation of the Hilbert transform is one of the crucial steps in our algorithm. This transform appears in computing the integrating factor in (8) and its Fourier modes in the angular variable (9). In order to perform a reliable numerical integration we need to properly account for the presence of the singularity in the kernel. The next simple lemma is key to our numerical treatment of the Hilbert transform.
where g, g a , g b are bounded and continuous functions defied by It means that the integrals on the right hand side is those in the sense of Riemann.
Proof. Firstly, for s ∈ (a, b), then since f ∈ C 1 the function g(s, t) is bounded continuous function on R 2 , and thus the calculation is followed by Finally, for s = a, then It is easily seen by calculation that both expressions are equivalent. Therefore we can theoretically choose any non-zero interval to evaluate the Hilbert transform as the Riemann integral of bounded and continuous functions. For the case of a < s < b, we split the interval at s in numerical computation: because g(s, t) may not be smooth on t = s. It is also convenient to employ the mid-point rule in order to avoid implementation of f appeared in g.

4.2.
Computation of Cauchy-type Integral Formula. Assume that ∂D has a parameterization ζ(ω), 0 ≤ ω < 2π and D contains the origin for simplicity. Suppose that 0 = ω 0 < ω 1 < · · · < ω K−1 < ω K = 2π. Take ζ k ∈ ∂D with Arg ζ k ∈ [ω k , ω k+1 ), where Arg ζ ∈ [0, 2π) is the argument of ζ ∈ C. Let us consider the discretization of the complex integral of an integrable function f (ζ) = We consider two examples as Ψ k . If Ψ k • ζ is the characteristic function of the interval [ω k , ω k+1 ], then the integral is approximated by where ζ k is the derivative at a certain point in the interval [ω k , ω k+1 ). On the other hand, if one can choose Ψ k • ζ be the continuous and piecewise-linear function with Ψ k • ζ(ω ) = δ k (Kronecker's delta) and ω k = Arg ζ k , then the trapezoidal rule gives an approximation as where ζ k = ζ (ω k ), ζ −1 = ζ K−1 − 2π and ζ K = ζ 0 + 2π by virtue of the periodicity. For the purpose of numerical computation, we express the integral in the second term on the right hand side of (6) as with F being the series in (6).
Step 1. Fix positive integers M and S. The integer S should be chosen sufficiently large, at least S ≥ M +3. Let ζ k = ζ (Arg ζ k ), 0 ≤ k < K. Choose the segmentation 0 = ω 0 < ω 1 < · · · < ω K = 2π so that Arg ζ k ∈ [ω k , ω k+1 ). From the periodicity we assert that ω j+K = ω j + 2π, j ∈ Z. We introduce an inscribed polygonal domain D h ≈ D, and take a triangulation T = { τ } of D h , i.e. each τ is a triangular domain, τ ∩ τ k = ∅ if = k, and D h = τ . Let P 1 (T ) denote the set of the piecewise linear continuous functions with respect to T . We denote by V the set of vertices of T .
Step 2. Compute and for 0 ≤ s ≤ S, 0 ≤ k < K, and z ∈ V . The function h[µ t ] is evaluated by the use of mid-point rule as stated so far.
Step 7. Step 8. For m = M + 1, M, . . . , 1, 0 (in descending order), find a piecewise linear continuous function I m | ∂D = a ϕ by interpolating { I m,k ; 0 ≤ k < K } obtained in Step 3, where ϕ is a periodic and piecewise linear continuous function on ∂D with ϕ (v k ) = δ k (Kronecker's delta), v k ∈ V ∩ ∂D. A more detailed example follows the algorithm description.
Step 9. For m = M −1, M −2, . . . , 1, 0 (in descending order), find an approximation I m ∈ P 1 (T ) to I m ∈ H 1 (D) by solving the Dirichlet problem of the Poisson equation (13) with the standard P 1 finite element method [11]. The variational formulation for (13) is approximated as follows; Find I m ∈ P 1 (T ) with (13b) so as to satisfy Firstly for m = M − 1, Step 7 and Step 8 give I M +1 and I M on the right hand side at z ∈ V , which leads the interpretation I M +1 , I M ∈ P 1 (T ). Particularly, if I M +1 (x) = a x 1 +b x 2 +c on a triangle τ , then ∂I M +1 | τ = (a −ib )/2. Similarly for the test function ϕ ∈ P 1 (T ), we can find ∂ϕ. The integration on the right hand side can be evaluated by a Gauss-type numerical integration [11] on each triangle τ . Then we can obtain I m−1 ∈ P 1 (T ) by solving the linear system.
Step 10. For each triangle τ ∈ T , let I 1 | τ = a x 1 +b x 2 +c . Then reconstruction of q| τ is given by (12a) as This ends the Algorithm. In Step 8, there may be a mismatch between the measurement points on the boundary and vertices of the triangulation for reconstruction. We solve this mismatch by interpolating the data on formers. Below we detail an example.
Let us assume that k I m,k χ k gives an interpolation of { I m,k ; 0 ≤ k < K }, where χ k is square integrable on ∂D. Then a 1 , . . . , a L can be determined as its best approximation in the sense of least square Hence the boundary value is given by I m (v ) = a . The similar strategy is applicable for Step 6.

4.4.
A locally optimal truncation criterion. We give a criteria on the choice of M . In the proposed algorithm, I 0 is obtained in complex values. On the other hand, the exact value of the zero-th Fourier mode I 0 is real-valued since I is so. Therefore the imaginary part of I 0 comes as errors in reconstruction. It is reasonable to consider that the errors in the real and the imaginary part interact each other. This observation leads a choice of M so as to minimize the imaginary part of I 0 , which is expected to reduce the error in the real part efficiently. Based on this consideration, we call M optimal which attains a local minimum of the imaginary part where |τ | is the area of τ and µ t (z ) − 2πµ s (z )p 0 (z ) Im I 0 (z ) corresponds to the imaginary part of reconstructed q on τ . In order to obtain an optimal M , we reconstruct the source for several values of M , then choose a value which minimizes E imag . We stress here that the optimality indicator in (18) does not require knowledge of the unknown source.

Numerical Experiments
In this section we demonstrate the numerical feasibility of the proposed algorithm for two numerical examples. All computations are processed with IEEE754 double precision arithmetic. In both numerical experiments, the measurement data is generated by solving the forward problem (1) by the piecewise-constant upwind approximation [13]. Therefore it is natural to use the characteristic function as the interpolation basis Ψ k in Section 4.2, and χ k in Step 8. The triangulation is generated by FreeFem++ [15].
The scattering coefficient is µ s ≡ 5. Physically, It means that the particle scatters on average every 1/5 unit of length. Given that D is the unit disc, particles scatter on average 10 times before getting out. We use the the two dimensional version of the Henyey-Greenstein scattering kernel in (2) p(ξ · ξ ) = 1 2π with the anisotropy parameter g = 1/2. This choice if half way between the ballistic g = 1 and isotropic g = 0 case. The Fourier expansion of the scattering kernel in (15) 1 2π yields the simple form of the modes p m = g m /2π, for all m ≥ 0.
In the discretization of the boundary, we adopt S = 128, while for the mid-point rule in the computations of the integral transforms, we use 100 sampling points.
be the rectangular, respectively circular subsets of D as shown in Figure 1. The source, to be reconstructed, is otherwise.
The absorption coefficient µ a (x) is given by To generate the boundary data, we solve the forward problem using the numerical method in [13] with a triangular mesh of 5, 542, 718 triangles, and 360 equispaced velocity intervals to describe the velocity directions. We disregard the value I(x, ξ) of the solution for x ∈ D, and only keep the boundary values. The obtained boundary data I(x, ξ) on ∂D × S 1 is depicted in Figure 2. In this figure, for x ∈ ∂D indicated by cross symbols (×), the graph of I(ζ, ξ), ξ , ξ ∈ S 1 are shown by a red (closed) curve in the polar coordinate with the center at each ζ. In other words, the red curve is the graph of { ζ + 2I(ζ, ξ)ξ ; ξ ∈ S 1 } for ζ ∈ ∂D indicated by the cross symbols. By the assumption of no incoming radiation, i.e. I| Γ− = 0, the red curve never appear inside |x| < 1 indicated by gray in this expression.
To better exhibit the effect of scattering, in the second representation of the data I Γ+ in Figure 3 we use the same coordinates as those used in a classical sinogram for the Radon transform data. More precisely, if ξ ⊥ ∈ S 1 denote the counterclockwise rotation of ξ by π/2, then the horizontal axis is the argument   Figure 2, the horizontal axis is the argument of the projection plane, and the vertical axis is the distance from the projected origin.
of the projection plane with direction ξ ⊥ , while the vertical axis is the distance from the projected origin. In the absence of scattering, p = 0, this representation would be exactly the sinogram for the attenuated Radon transform data. Note in Figure 3 how the scattering had combined and smeared out the features from the three locations of the source.
To avoid an inverse crime, the triangulation used in the reconstruction is different from that in the forward problem. In particular the reconstruction mesh consists of 6, 998 triangles with 5, 400 vertices (much less than the 5, 542, 718 triangles used in the forward problem), and is generated without any information of the location of the subsets R, B 1 , and B 2 . The computational time for reconstruction with M = 6 is 340 seconds on Xeon E5-2650 v4 (2.2GHz, 12 cores) with OpenMP. Almost all computational time are occupied by the computation of the discrete Fourier transform of e −h and e h in Step 2, and the boundary integral (17) in Step 5, as 117 seconds, 79 seconds, and 141 seconds respectively.
The reconstructed q(x) with M = 6 is shown in Figure 4. Its cross sections along the dotted diameters x 2 = − √ 3x 1 and x 2 = 0, passing through the origin and the center of B 2 and R respectively, are depicted in Figure 5. The reconstructed q(x) shows a quantitative agreement with the exact source in (19). Similar to the Xray and attenuated X-ray tomography, the artifacts appear due to the co-normal singularities in the source but also in the attenuation.
where |τ | is the area of τ , z is the center of the triangle τ , and the summation runs over the triangles where q is continuous (in particular, constant in the example). From the figures, the error in imaginary part scaled on the right axis is sufficiently smaller the pseudo-error scaled on the left axis. Both errors take minimum around M = 5, and increase after that. More precisely, the error of the imaginary part (18) is minimum at M = 6, while the pseudo-error of reconstructed q is minimum at M = 3. According to the criteria stated before, we adopt M = 6, which causes |p M | < 2.49 × 10 −3 in the setting. x along x 2 =0 exact Figure 5. Section of the reconstructed source along the dotted lines in Figure 4, |x 2 + √ 3x 1 | < 0.05 (on the left) and |x 2 | < 0.05 (on the right) respectively. The arrow on the right figure indicate ∂B 1 where the scattering coefficient is discontinuous.

Experiment 2.
Let us consider the situation that µ a is given by the modified Shepp-Logan phantom [36], which occupies the ellipse We consider three disjoint domains in D (Figure 7): and the internal source is given by otherwise. The boundary measurement is generated by solving the forward problem with 1, 554, 282 triangles and 360 velocity directions. The numbers of measurement points are K = 3, 000 on ∂D, and N = 360 on S 1 . On the contrary, the reconstruction mesh consists of 6, 010 triangles with 3, 106 vertices. Similarly as the previous experiment, the latter mesh is generated without any information of µ a and q. The boundary nodes on the ellipse generated by FreeFem++ are not equi-spaced with respect to their angles. Figure 9 depicts the reconstructed q(x) on D with M = 8, while Figure 10 is its sections on the dotted lines. The computational time for reconstruction with M = 8 is 247 seconds on Xeon E5-2650 v4 (2.2GHz, 12 cores) with OpenMP. From the results, the support of q is clearly and quantitatively reconstructed, while the profile of µ a does not appear.
The pseudo-errors in the reconstructed source and the errors of the corresponding imaginary part are shown in Figure 11. Our proposed optimality criterion yields to choose M = 8 and M = 40 as reasonable orders of truncation. Figure 12 shows sections of the reconstructed q(x) with M = 40 on the dotted lines in Figure 9. Although the value of (18) at M = 40 is smaller than that at  Figure 9, |x 2 | < 0.01 (on left) and |x 2 + 0.605| < 0.01 (on right). M = 8, the peak of the source in Q 2 is not obtained well with M = 40. The possible reason for this fact is that the size of Q 2 is small relative to the size of the triangular mesh, yielding that the discrete L 2 -norm used in computing the pseudo-error in (20) be less effective.
In general, the choice of an optimal truncation parameter M is not clear. However, our algorithm includes an optimality criterion which is independent of the knowledge of the source, thus making it feasible. The numerical experiments based on this choice were shown to produce accurate reconstructions.