Forward and Inverse Problem Formulation of Optical Tomography Based on Equation of Radiative Transfer

Optical tomography is non-invasive diagnostic technique. Mathematically, it is related to the evaluation of optical parameters from the equation of radiative transfer with diffused boundary measurements. Since the radiative transfer equation in expressed in the form of phase space, it is quite challenging to solve it computationally. In this paper, reconstruction is based on equation of radiative transfer in frequency domain and termination criteria for forward problem is proposed as ratio of residues. Directional and spatial variables of equation of radiative transfer are discretized with discrete ordinate method and ﬁnite volume method respectively. The sparse structure of matrices of complex valued algebraic equations is formulated.


Introduction
Optical tomography is a non-invasive medical imaging method utilizing near infrared light to explore biological tissue in order to estimate qualitative or quantitative information about the optical properties of the tissue. These optical properties can be employed for diagnostic purposes. Applications of optical tomography include, for example, brain imaging [1] [2] [3], breast imaging [4] [5] [6] and imaging of joints [7] [8] [9].
Mathematically, it is related to the computation of optical parameters for the equation of radiative transfer.
The confined validity of diffusion approximation only to specific contexts is proved by numerous theoretical and experimental analyses. For example, when the scattering features are more significant than the absorbing features [10].
Hence, there are numerous applications for which equation of radiative transfer based reconstructions are required. Examples are brain imaging, where the low absorbing and low scattering cerebrospinal fluid corresponds to non-diffusive light propagation; imaging of joints [9], where clear synovial fluid possesses non-diffusive light propagation [11], and small animal imaging [12], where small object dimensions (1-2 cm in diameter) elaborate the drawbacks of the orthodox diffusion approximation.
Algorithms are devised and experimentally proved in biomedical imaging for two and three dimensional problems employing a time independent and time dependent equation of radiative transfer [13] [14] [15] [16]. Though these innovations are a vital achievement towards practical implementation, but still appreciable cross talks are observed between optical reconstructions.
The cross talk means that purely scattering features are mainly mapped with absorption features or vice versa.
To overcome cross talks, which may result in wrong diagnosis, variety of data is required [17]. This is made possible in the recent years by the application of frequency domain measurements. Since frequency domain measurements provide information about both the phase and intensity of the waves, so better separation of scattering and absorption effects is expected [17] [18].
The paper is organized as follow. In section II, the problem of optical tomography is formulated in terms of equation of radiative transfer in frequency domain. In section III, discretization of equation of radiative transfer is presented. Sparsification and Modified GMRES are devised in section IV. Adjoint method for gradient computation is described in detail in section V. The implementation of the minimized procedure based on quasi -Newton algorithm are covered in section VI. Concluding remarks are offered in section VII.

Problem Formulation
The equation of radiative transfer in frequency domain describes the density of photon in the phase space. The phase space is a function of position x ∈ D ⊂ R n and direction θ ∈ S n−1 (Unit sphere in R n ) and equation is given by is unknown quantity and is the radiant power per unit solid angle per unit area perpendicular to direction of propagation at x in the direction θ. dµ is the surface measure on S n−1 . k(θ, θ ) is the scattering kernel that describes the probability of photons travelling in direction θ scatter into direction θ. It is a positive function independent of x and satisfies normalization condition. It is chosen as Henyey-Greenstein phase function due to highly peak forwardness of tissues given in [19]. In optical tomography, it is assumed that no photon travels in an inward direction at the boundary ∂D and boundary sets are given in [19].
Restrictions are imposed for optical parameters, i.e., σa and σs in order to avoid negative values by the introduction of space parameter given in [19].
As a result, the forward problem is well posed and a unique solution u(x, θ) is attained. Limited information is extracted from the mapping of the incoming flux at the boundary into the outgoing flux because only angular averages of outgoing flux are attainable conventionally. Mathematically, it means that the spatial resolution is highly confined. This is the main hurdle behind the employment of uniqueness and stability constraints of inverse transport theory [20]. As consistency is required with the accessible evaluation technologies so an evaluation operator with slight modification is acquainted and is defined in [19] [21].
The estimation of the unknown absorption and scattering coefficients (σa,σs) ∈ Q within the domain D when the distribution of the light sources and the evaluated data z on ∂D are given is known as the inverse problem and is given by Generally, the inverse problem is severely illposed. To overcome ill-posedness, uniqueness of reconstruction should hold. Solution of (2.2) can be obtained by solving following least square formulation.
Extra evenness limitations are imposed on the coefficients to stabilize the least square problem by the introduction of space parameter Q ad given in [19]. The following regularized least square functional can be introduced [21].
where the last term is called a regularization term and α is the regularization parameter. The L-curve method is applied to choose the optimal regularization parameter α in Fα(σa, σs) [19] [22] [23]. Search for a pair (σa, σs) that minimizes the least square functional Fα(σa, σs). Fα(σa, σs) should have atleast one minimum. However the uniqueness of reconstruction is not guaranteed because Fα(σa, σs) is not strictly convex.
Inverse problem in optical tomography is implemented with a minimization approach relied on gradient evaluation. It is, therefore, required to evaluate the Fréchet derivative of the least square functional Fα(σa,σs) [19] [21]. However, it is not economical to evaluate the Fréchet derivatives directly because the optical features are infinite dimensional. Therefore, the adjoint method is adopted for the differentiation.

Discrete Ordinate Method
A discrete ordinate method was developed by Chandrasekhar in 1950 [24]. In this method [25], total scalar flux is approximated and is described as the integral of u(x, θ) over the S n−1 by following quadrature rule.
where θj is the jth direction and nj is the associated weight.

Finite Volume Method
Cell centred version of finite volume method is implemented to discretize spatial variables [19] [21]. Consider a mesh of R n , M , consisting of polyhedral bounded convex subsets of R n . This mesh M covers the computational domain D. Let C ∈ M be a control cell (an element of mesh M ). ∂C is its boundary. Vc is its Lebesgue Measure. Now the unknown quantity u(x, θj) is assumed constant in control cell C and is defined as u(x, θj) ≡ u c j on C. Integration over cell C along with the application of Divergence theorem, Lebesgue surface measure on boundary of cell C and result in the following expression.
where nC (x) denotes the outward normal to ∂C at point x ∈ ∂C, dγ(x) denotes the surface Lebesgue measure on ∂C. The flux S C,i θj · nC (x)ujdγ(x) can be approximated by Firstorder upwind technique given in [19] [21]. Hence the full discretization of discrete ordinate equation is given as:

Sparsification
The discrete transport equations are collected on all control cells. As a result, following system of complex-valued algebraic equations are obtained.
where X ∈ C N J×N J is discretized streaming collision operator, Y ∈ C N J×N J is discretized streaming scattering operator, Z is discretized boundary source f (x, θ), U ∈ C N J×1 consists of the values of u(x, θ) on the cell C in the direction θj and is systematized as: U = (u1, ..., uJ ) T with U j = u 1 j , ..., u N j T ∈ C N . The matrices X and Y are sparse block matrices. X is a block diagonal matrix that can be written as: where the discretization of the advection operator X is represented by Xj ∈ C N ×N and is defined by X u := θj · ∇u. From a first order upwind scheme, Aj does not contain more than N × NE non zero elements. NE denotes the total number of edges (surfaces in three dimensions) of each control cell. Matrix W0 ∈ C N ×N is diagonal: The matrix Y can be expressed as direct product of two smaller matrices: with E0 ∈ C N ×N a diagonal matrix given by, and K ∈ C J×J a dense matrix with component In realistic applications, the number of orientations J the number of constituents of spatial mesh N . The matrix K is not symmetric unless ηj is assumed to be constant. Hence, the matrix X − Y is neither symmetric nor positive definite.

Modified GMRES Algorithm
A preconditioned modified generalized minimal residual algorithm algorithm is chosen to solve the forward problem [26] [27] [28]. The modified GMRES algorithm is terminated if the relative residues are sufficiently small. The termination criteria is where U0 is the initial guess, U k is the U value at the k-th GMRES iteration. Termination criteria based on the ratio of residues increases the accuracy of data obtained from modelled parameters in forward problem. Relationship of guess measurements and measurements obtained after performing k iterations is determined as this ratio sets the criteria for true values.

Discrete Adjoint Problem
The discrete adjoint problem is designed to compute the gradient of discrete objective function with respect to the optical properties on each cell [21]. For simplification, let σa ∈ R N ×1 denotes the absorption coefficient vector, i.e., (σ 1 a , ..., σ C a , ..., σ N a ) T . Let σs ∈ R N ×1 denotes the scattering coefficient vector, i.e., (σ 1 s , ..., σ C s , ..., σ N s ) T . The discretized objective function required to minimize is of the form given as; (5.1) where number of detectors employed for each source are N d , z δ d denotes the calculation of the source, δ denotes the level of noise in the calulations. P d ∈ R 1×N is evaluation operator in discrete form that averages outgoing flux over S n−1 + . The regularization term in discrete form is given in [19]. Differentiating (5.1) with respect to simultaneously, differentiating (4.1) with respect to σ C a and then rearranging results in Now a new state variable V ∈ C N ×1 called adjoint variable of U and adjoint solution of (4.1)is given by Rearranging (5.5) and substituting in results in Similarly, differentiating (5.1) with respect to σ C s yields respectively. This is because, these matrices have very simple sparse structures according to (4.2) and (4.4). Instead of this, a matrix free method is implemented given in [19] [21].
Hence, ∂Fα ∂σ C a and ∂Fα ∂σ C s can be evaluated without construction of any intermediate matrices.

Numerical Implementation
The Quasi Newton optimization algorithm is implemented to compute the regularized least square problem (RLS). In realistic applications, the convergence rate is much faster than non linear conjugate gradient method with the employment of The Fletcher-Reeves or the Polak-Ribiére [14] [29]. Gauss Newton method without the confinement on constraints generally fails in frequency domain reconstruction because the extreme non linearity slows the convergence rate [29] [30].   With σ denoting the vector of discretized optical properties, the quasi Newton methods can be characterized by the following iterative process.
σ k+1 = σ k + α k p k (6.1) where k ∈ N. p k is a descent direction. The BFGS algorithm chooses p k to be the solution of an approximated solution of Newton type optimality equation p k = H k g k , α k is a step length, g k is the gradient of the least square functional and is given by; g k = −∇σFα(σ k ), H k is the inverse Hessian matrix of Fα at step k. The computation of real inverse Hessian matrices is very time consuming. Hence, the limited memory BFGS algorithm chooses to approximate H k by the following updating rule [29].
where W k = I − ρ k y k s T k and I is the identity matrix, s k = σ k+1 − σ k , y k = g k+1 − g k , ρ k = 1 y T k s k . To enforce limitations on optical parameters, the relation p k = H k g k is changed little by the implementation of a gradient projection method [29] [31] [32].

Conclusion
A better optical image reconstruction method is devised. This method considers all the possible errors related to the collection of data. Computation time is significantly reduced by devising sparse structures of complex valued algebraic matrices as low significant data is ruled out. Termination criteria based on the ratio of residues increases the accuracy of data obtained from modelled parameters in forward problem.The frequency domain reconstruction reduces crosstalk between absorption and scattering parameters hence the shortcomings of the diffusion approximation of the radiative transfer equation are overcome. It provides a very useful diagnostic approach in modelling strongly absorbing regions such as large blood filled spaces e.g brain haematoma, low scattering void like inclusions such as spaces filled with cerebrospinal fluid, amniotic fluid or synovial fluid and optically relatively thin media such as fingers and small animals imaging.