Solution of Radiative Transfer Equation with a Continuous and Stochastic Varying Refractive Index by Legendre Transform Method

The present paper gives a new computational framework within which radiative transfer in a varying refractive index biological tissue can be studied. In our previous works, Legendre transform was used as an innovative view to handle the angular derivative terms in the case of uniform refractive index spherical medium. In biomedical optics, our analysis can be considered as a forward problem solution in a diffuse optical tomography imaging scheme. We consider a rectangular biological tissue-like domain with spatially varying refractive index submitted to a near infrared continuous light source. Interaction of radiation with the biological material into the medium is handled by a radiative transfer model. In the studied situation, the model displays two angular redistribution terms that are treated with Legendre integral transform. The model is used to study a possible detection of abnormalities in a general biological tissue. The effect of the embedded nonhomogeneous objects on the transmitted signal is studied. Particularly, detection of targets of localized heterogeneous inclusions within the tissue is discussed. Results show that models accounting for variation of refractive index can yield useful predictions about the target and the location of abnormal inclusions within the tissue.


Introduction
A special attention in diffuse optical tomography is focused on the development of methods for detection of photons providing the information concerning optical parameters of the explored medium. This gives the targets of localized nonhomogeneous inclusion arising in tissues due to various pathologies, like tumor formation, local increase in blood volume, and other abnormalities [1][2][3][4]. In radiative transfer theory, the most used parameters in modeling laser radiation interaction with biological tissue are absorption and scattering [5][6][7]. However some other studies evoked a significant variation of refractive index of abnormal biological tissues especially in the near infrared range. More precisely, experimental results [8,9] showed that the tissue of malignant tumors could manifest an increase of the refractive index which can attain until 10% of that of a normal tissue which encircles them. So, medical imaging by diffuse optical tomography should take advantage from the emergence of a third contrast parameter which is the refractive index. This led to the appearance of a big number of numerical and fundamental works in the field of radiative transfer in a varying refractive index biological medium. While the conventional radiative transfer equation (RTE) has been widely used to study interaction of near infrared radiation with biological media, there exist a number of works dealing with a modified radiative transfer equation in spatially varying refractive index media [10,11]. Some of these papers are interested in varying refractive index biological tissues [12][13][14]. In the present paper, our first concern is to contribute to the usability of the radiative transfer theory in a potential optical tomography setting in medical imaging. At this level, studying the effect of refractive index on the transmitted light through a biological rectangular layer should be crucial. This could improve detectability of heterogeneous objects in a typical tomography scheme. However, it is important to note that in a varying refractive index medium, the rays are not straight lines but curves. So even in a rectangular geometry, the varying index radiative transfer equation displays the classical form of the angular derivative terms commonly appearing when dealing with spherical and cylindrical geometries with uniform refractive index [15][16][17]. This finding gives rise to the use of Legendre transform as a manner for modeling these terms. Although this technique was used by Sghaier et al. [17] in a uniform refractive index spherical domain as an innovative view to handle these terms, it prevails as useful in this contemporary problem. This fact is our second concern in this paper. So, we present a computational RTE-based model suitable for basic diffuse optical tomography forward problem with spatially varying refractive index biological medium. We treat angular derivative terms by using the Legendre integral transform technique. We investigate cases concerning optical tomography applications. Results concerning the effect of the refractive index variation on the detected signal are shown.

Mathematical Model
In this work, the radiative transfer equation in a human biological tissue is described by using a stationary varying refractive index RTE [18,19] where (i) Ψ( ⃗ , ⃗ Ω) is the directional energetic radiance at the spatial position vector ⃗ = ( , , ), (iii) ( ⃗ ) and ( ⃗ ) are the absorption and scattering coefficients, respectively, (iv) = vac / is the ratio of speed of light in a vacuum, vac , and the refractive index , (v) the source term ( ⃗ , ⃗ Ω) is an injected radiance at the medium's boundary, (vi) the phase function ( ⃗ Ω, ⃗ Ω ) describes the probability that, during a scattering event, a photon with direction ⃗ Ω is scattered in the direction ⃗ Ω.
Equation (1) takes into account the fact that the rays are not straight lines but curves. It involves terms that illustrate the expansion or the contraction of the cross section of the tube of light rays in the medium.
For a two-dimensional problem and in Cartesian coordinate system of the -plane, the terms due to the refractive index variation can be expressed as where cos and sin are the Cartesian coordinates of the unit direction vector in the -plane. In fact we assume that the radiance of out of plane directions is negligible. By using notations = cos and = sin , (2) displays the classical form of the angular redistribution term commonly appearing when dealing with spherical and cylindrical geometries with uniform refractive index [20] The angular redistribution terms will be noted so (3) becomes

Numerical Method
In our numerical implementation, we use a rectangular domain which is divided into a set of × elementary uniform volumes Δ = Δ Δ Δ with a uniform unitary depth (Δ = 1). The angular discretization is obtained through a discrete ordinate technique. This yields a set of discrete directions, , = 1 ⋅ ⋅ ⋅ giving a set of angular discrete direction cosines ( , ), = 1 ⋅ ⋅ ⋅ . An orientation depending on the incident ray direction is adopted for each cell [7]. Calculations are done by using integration of (1) over an elementary volume Δ for each discrete direction. This gives Computational and Mathematical Methods in Medicine   3 where , and , are the discrete angular derivative terms at the angular ordinates and , respectively, and is a weighting factor. The discrete term of Henyey-Greenstein phase function is written as where is the anisotropy factor. If the direction cosines are positive, the directional radiance is known on the faces and and they are unknown on the faces and of the ( , )-cell and also in the centre . Therefore, we need two complementary relations to eliminate Ψ , and Ψ , ; this can be obtained by using interpolation formula where is an interpolation parameter. Using these relations, becomes Theoretically, if we know the solution in the ( , )-cell, we can do calculus over the cells ( + 1, ) and ( , + 1) using the boundary conditions and the following relations: If the direction cosines are both positive, we get the following equation: this gives

Numerical Treatment of Angular Derivative Terms with
Finite Legendre Transform. As is explained in [17], we consider the following Legendre transforms: where is the degree Legendre polynomial. According to Sturn-Liouville equation defining Legendre Polynomials, we can write To obtain derivative terms, we make use of numerical quadrature: The above systems of equations are closed by the obvious relations: So discrete derivative terms in the ( , )-cell into the medium can be obtained by solving the linear algebraic equationŝ where A = (  . . . and ,̂, , , , , are straightforward by replacing by . A set of weights are judiciously chosen with an equally spaced distribution of angular ordinates ( , ). The matrices and must be regular and they are inversed each only once in all the calculation process.
To solve (12), we use successive iterations to actualise the implicit internal source term in the right member. So, this gives The iteration process is repeated until a convergence criterion is attempted. To improve convergence speed, we use a successive overrelaxation method. So the updated value Ψ +1 , , is a linear combination of the iterated value Ψ , , and the previously computed value where is a relaxation parameter whose value is usually between 1 and 2. The solution is obtained when the relative discrepancy value is smaller than a tolerance value. In all our calculus, we have taken 10 −8 as a tolerance value. As initial condition, we take a field of null radiance. Also, all our calculations are done in the case of interpolation diamond scheme ( = 0.5).
If the direction cosines are not both positive, the precedent equations are valid provided that the orientation WESN of cells is done according to the direction of propagation [7]. In all our investigations, the injected power source is assumed to be equivalent to a forward collimated monochromatic intensity placed at a source point on the middle of the bottom side of the boundary. Results shown below are obtained by using a continuous wave source with a uniform equivalent intensity value of 50 mW⋅cm −1 .

Boundary Conditions.
On the boundary, the radiance is the sum of the external source contribution and the partly reflected radiance due to the refractive index mismatch at the boundary where ⃗ is a position on the boundary and ⃗ is an outer normal unit vector. The reflectivity can be calculated for each direction using Fresnel's relations. Medicine   5 To present our results, we use the detected fluence rate which is given in a ( , )-detector point on the boundary as

Computational and Mathematical Methods in
with , , Also, we make use of a normalized detected fluence rate defined as where is the number of the detector points on one side of the boundary and is a weighting factor from the generalized trapezoidal integration rule.
In all calculations, we have used 28 detector points on each side. Also, all calculus is carried out by using 16 uniformly distributed discrete directions and a space grid of 121 × 121 cells.

Continuous Varying Refractive Index Medium.
In this investigation, we study near infrared radiation transport in a rectangular medium exposed to a continuous collimated source which is placed on the bottom side of the boundary. Figure 1 shows the considered medium; it is assumed to be 2 × 2 cm sized with varying refractive index. Within the medium, we consider a -axis linear refractive index variation with different gradient values. To show the effect of the refractive index on detected fluence rate, we have used a weakly absorbing and forward scattering background medium whose optical parameters are shown in Figure 1. Figures 2 and 3 represent the response of the medium through the detected signal on the top and right side of the boundary in the case of linear variation of refractive index. We note that the response of the medium varies linearly according to the gradient of refractive index. The maximum transmission moves to regions of height index ( Figure 2). Indeed, according to boundary conditions of the medium (Descartes Laws) the transmission window increases proportionally of refractive index between the medium and outside medium (air). On the other hand, the right side ( Figure 3) shows the response of the medium detected on the right side of the boundary in the case of linear and parabolic variation, respectively. The detected fluence rate curves presents distinguished effect of refractive index gradient in linear case. The transmission zone on the boundary increases with decreasing gradient values. Even though most detected transmission is obtained near the source, a weak gradient can augment transmitted radiation relatively far from the source while high values of gradient can block transmitted radiation within the medium.

Effect of Stochastic Varying Refractive Index.
In this investigation, we keep the same medium precedent, but we will study the stochastic variation refractive index. We control   this variation by uniform random fixing of a given interval. The other optical properties are the same as in the precedent investigation. In such cases, there is an obvious effect of the heterogeneity on the detected signal especially when the refractive index is increased or decreased by 10% in a stochastic way. The detected signal on different sides of the medium is shown in Figures 4 and 5. There is a significant effect of this variation on the detected signal on the top side and right side, because there is strong perturbation in medium. These findings highlight the potential of refractive index as a possible detection parameter of a tumor in a surrounding safe tissue. These findings open prospects for further study

Conclusion
This study attempted to develop a computational way helping in detection of abnormalities in a biological tissue. This should enable predictions of eventual tumor existence when using a diffuse optical tomography scheme. The used model is based on stationary radiative transfer equation including a possible continuous and stochastic variation of refractive index. In particular, computational technique of Legendre transform is extended to handle angular derivative terms arising by the varying refractive index consideration. The obtained computational model is implemented to investigate some practical situations in diffuse optical tomography (DOT) setting. Obtained results showed that variation of refractive index can yield useful predictions about the target and the location of abnormal inclusions within the tissue. These findings open prospects for further study of the effect of local refractive variation in time-domain and frequencydomain schemes in diffuse optical tomography.