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.

1. 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-4]. In radiative transfer theory, the most used parameters in modeling laser radiation interaction with biological tissue are absorption and scattering [5-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-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-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.

2. 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

is the directional energetic radiance at the spatial position vector ,

is the refractive index distribution,

and are the absorption and scattering coefficients, respectively,

c = c vac/n is the ratio of speed of light in a vacuum, c vac, and the refractive index n ,

the source term is an injected radiance at the medium's boundary,

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 x-y 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 x-y 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

3. Numerical Method

In our numerical implementation, we use a rectangular domain which is divided into a set of I × J elementary uniform volumes ΔV = ΔxΔyΔz with a uniform unitary depth (Δz = 1). The angular discretization is obtained through a discrete ordinate technique. This yields a set of M discrete directions, φ , m = 1 ⋯ M giving a set of angular discrete direction cosines (ξ , η ), m = 1 ⋯ M. 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 ΔV for each discrete direction. This gives where D and D are the discrete angular derivative terms at the angular ordinates ξ and η , respectively, and w is a weighting factor. The discrete term of Henyey-Greenstein phase function is written as where g is the anisotropy factor. If the direction cosines are positive, the directional radiance is known on the faces W and S and they are unknown on the faces E and N of the (i, j)-cell and also in the centre P. 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, (6) becomes

Theoretically, if we know the solution in the (i, j)-cell, we can do calculus over the cells (i + 1, j) and (i, j + 1) using the boundary conditions and the following relations:

If the direction cosines are both positive, we get the following equation: this gives

3.1. Numerical Treatment of Angular Derivative Terms with Finite Legendre Transform

As is explained in [17], we consider the following Legendre transforms: where P is the m 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 (i, j)-cell into the medium can be obtained by solving the linear algebraic equations where and are straightforward by replacing ξ by η.

A set of weights are judiciously chosen with an equally spaced distribution of angular ordinates (ξ, η). The matrices A and A 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 Ψ 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.

3.2. 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 boundarywhere is a position on the boundary and is an outer normal unit vector. The reflectivity R can be calculated for each direction using Fresnel's relations.

;

and ,

To present our results, we use the detected fluence rate which is given in a (i , j )-detector point on the boundary as with

Also, we make use of a normalized detected fluence rate defined as where D is the number of the detector points on one side of the boundary and w ′ 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.

4. Results and Discussion

4.1. 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 x-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.

Figure 1

A test-medium with background properties and detector points.

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).

Figure 2

Response of a continuous varying refractive index medium detected fluence rate on the top side: gradient effect.

Figure 3

Response of a continuous varying refractive index medium detected fluence rate on the right side: gradient effect.

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.

4.2. 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 of the effect of local refractive variation in domain and frequency-domain schemes in diffuse optical tomography.

Figure 4

Response of stochastic varying refractive index medium on the top side.

Figure 5

Response of stochastic varying refractive index medium on the right side.

5. 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 frequency-domain schemes in diffuse optical tomography.