1. Introduction

Optical tomography uses light for visualization of three-dimensional structure of transparent and semi-transparent objects. An application of the inverse Radon transform to optical projection data allows to obtain cross-sectional high-resolution images of small non-scattering objects such as an embryo or foetus. Optical Projection Tomography (OPT) was developed as an optical analogue of X-ray Computerized Tomography (CT) [1, 2] and has found many important applications as an imaging methodology in medicine and biology [3–6]. OPT is likely to evolve to a powerful and sophisticated imaging modality by taking advantage of time-gated technology [7], structured illumination patterns, and light polarization. On the other hand, in contrast to CT, it requires correction for unknown refractive index inside an object. Presence of scattered light in a sample introduces even more complexity to this technique.

Light scattering can be addressed by using Monte-Carlo methods [8–10] or employing the Radiative Transfer Equation (RTE) [11–17]. Both methodologies are computationally very expensive. Image reconstruction algorithms involve repeated solution of the direct problem and, therefore, the use of Monte-Carlo or the RTE significantly affects the performance of these methods. Recently it was suggested to account for only singly scattered photons with angularly selective intensity measurements and apply the broken ray Radon transform [18, 19]. This approach is applicable to weakly scattering media but is invalid in the presence of highly scattering inclusions.

In this paper we consider an approximation to the RTE, which is based on the assumption that scattering media consist of weakly and highly scattering regions, whose transport coefficients differ by an order of magnitude. Such a situation is quite common in Nature. For instance, the atmosphere can be considered mostly non or weakly scattering with clouds presenting highly scattering inclusions. A live embryo or foetus is mostly transparent with internal organs, brain or eye balls being highly scattering.We also make an assumption that the phase function, which is included in to the RTE collision term, can be approximated by the first two terms in its Fourier series expansion. This assumption lets us use the telegraph equation approximation for the light transport inside highly scattering inclusions. For the present, we neglect refractive index variation in the medium.

The image reconstruction approach considered in this paper employs angularly selective intensity measurements. Thus, optical parameters of a scattering medium are reconstructed by using optical projection datasets and scattered light outgoing from the medium at some angles with respect to the direct radiation. The direct radiation is assumed to take the form of parallel rays entering the medium. The image reconstruction algorithm is based on the variational framework and involves repeated numerical solution of first and second order differential equations in the Fourier domain.

This paper is organized as follows. In the next section, we discuss the direct problem and derive an approximation to the RTE, which is valid under our assumption. Implementation details of the direct solver and numerical simulations illustrating our approach are presented as well. At the end of this section the reconstruction algorithm is described. Section 3 is devoted to numerical experiments where reconstructions are presented and discussed.

2. Methodology

2.1. Direct problem

Light transport in scattering media is modeled by the RTE, which is the integro-differential equation describing a balance of radiation along an given direction s [20]. In this study we treat the RTE as a first order partial differential equation

where the function B in the source term reads

In Eqs. (1) and (2), (i) I denotes the intensity of light; (ii) μ̃ = μ + iω/c is the complex extinction coefficient, where μ is the transport coefficient, which is a sum of the scattering, μ_s, and absorption, μ_a, coefficients, ω is the Fourier parameter, and c is the speed of light; (iii) λ denotes the albedo of a single scattering event such that λ ∈ [0,1], μ_s = λμ and μ_a = (1−λ)μ; (iv) p(s · s′) is the phase function. The simplest form of an anisotropic phase function is assumed in this study

where ε ∈ [−1,1]. The last term in Eq. (2) is scattered once direct radiation I ₀ (r, s ₀). Note that the transport coefficient, μ, and albedo, λ, are chosen instead of conventional scattering and absorption coefficients. The albedo is a photon's probability to survive a single scattering event and, therefore, it controls the true absorption. The transport coefficient, which is reciprocal to photon's mean path length between successive scattering events, describes scattering properties of the medium. When the scattering medium is a mixture of two types of particles with different scattering properties then the resulting transport coefficient and the albedo are found as μ = μ ₁ + μ ₂ and λ = 2λ ₁ λ ₂ (λ ₁ + λ ₂)⁻¹, where μ_j and λ_j (j = 1, 2) are corresponding parameters for each type of particle. Parameters of a mixture of more than two types of scattering particles are computed recursively.

The direct radiation serves as a source for scattered light and is assumed in the form of parallel rays entering the domain along the unit vector s ₀. The intensity of each ray (the Green function) is found by solving the equation

where Q ₀ is the source amplitude and r ₀ belongs to a source plane r ₀ · s ₀ = const. The solution of Eq. (4) is given by

Sometimes it is more convenient to integrate from the observation point r rather than from r ₀. The source and observation points are connected by a line r = r ₀ +s ₀ l, and, therefore, the integration in exponent in Eq. (5) can be rewritten in the form

∫^l ₀ μ̃(r ₀+s ₀ l′)dl′ = ∫^l ₀ μ̃ (r−s ₀ l′)dl′.

Finally, the general solution of the RTE at the observation point r is given by

The direct radiation I ₀, Eq. (5), must be added to I when s = s ₀. To correspond to physical reality we add the direct radiation to the scattered intensity when s · s ₀ ≥ 1−δ, where δ is finite but small number.

The intensity, Eq. (6), can be computed numerically when the function B is known. Exact knowledge of the function B is equivalent to computing the solution of the RTE. Here we suggest an approximation to the function B according to the assumption that the medium consists of weakly and highly scattering regions, whose transport coefficients differ by an order of magnitude. We further assume that recorded photons coming from weakly scattering regions are scattered only once, i.e. 1/μ is a length scale on the order of physical dimensions of the scattering domain. This assumption implies a presence of photons scattered more than once. However, they do not reach the CCD array. Therefore, the method of successive approximations applies [20] when only the first approximation is retained. Thus, as an approximation for B we take the scattered once direct radiation p(s · s ₀) I ₀ (r, s ₀), which is the last term in Eq. (2). On the other hand, in highly scattering regions the intensity in Eq. (2) is approximated by

where u denotes the average intensity defined by

and κ is the diffusion coefficient

The expression for the approximate intensity, Eq. (7), results from first two moments of the RTE. That is, consequent multiplication of Eq. (1) by 1 and by s and integration over the whole solid angle leads to a system of two first order differential equations for u and the flux q = −κ∇u. This system is closed when the phase function is assumed in the form of Eq. (3). Elimination of the flux in this system results in the Helmholtz equation for the integrated intensity [21]

where

and the source term, ρ, represents direct intensity averaged over the whole solid angle [21]. Substitution of Eq. (7) into Eq. (2) simplifies the function B to

This final form of the term B is used for computing the intensity I, Eq. (6), where in weakly scattering regions the integrated intensity, u, is neglected.

2.2. Implementation details

High resolution imaging imposes certain constraints on mesh density. A dense mesh requires high performance algorithms for solving the direct and inverse problems. In many cases high performance is achieved by using the simplest and computationally inexpensive approaches. In order to spare computational resources, efficient dynamic memory allocation is the one of the high priority issues. In this study a Cartesian mesh has been chosen. The entire computational domain is split into computational cells (voxels), whose dimensions correspond to a pixel's dimension of the CCD array. All functions are approximated by piecewise constant functions having constant values in each computational cell. The Helmholtz equation (Eq. (10)) is solved by applying the Finite Volume method, which is identical to the Finite Difference numerical scheme for the regular Cartesian grid. However, this approach can be further generalized for adaptive structured meshes.

Computation of the intensity in the scattering medium involves integrations along the rays' paths. A ray integration in Eq. (6) is performed by using Siddon's algorithm [23], which is the ray tracing algorithm designed for Cartesian grids. It finds distances between successive intersections of a ray with coordinates planes. These distances, together with cell's indices, are computed for planes parallel to the x, y, and z-axes. The final three-dimensional ray's path is found by applying a merge sort to these three sets of distances.

The choice of anisotropic phase function in the form (Eq. (3)) requires numerical integration of ελμκ s·∇u along a ray in accordance with Eq. (12). For the sake of simplicity a constant value of the parameter ε is assumed everywhere in the domain. Then, we are looking for an inexpensive way of numerical evaluation of the following term

Firstly, we note that a piecewise constant function u satisfies

where ${\left[u\right]}_{l_i}$ denotes a jump of u across a cell's interface at l = l_i along the direction s. Then, the line integral (Eq. (13)) is approximated by a sum

where ${\lambda }_{l_i}$ and ${\mu }_{l_i}$ are the cell's interface values of the albedo and the transport coefficient at l = l_i. The diffusion coefficient at the cell's interface, ${\kappa }_{l_i}$ , is computed according to Eq. (9) where the interface values of λ and μ are substituted. The distance Δl_j is the length of the ray's path within a cell provided by Siddon's algorithm. Here, the index j enumerates cells on the ray path and μ̃ is the extinction coefficient of j-th cell. The cell's interface values of parameters ${\lambda }_{l_i}$ and ${\mu }_{l_i}$ are set to the averages

{λ} = (1/2) (λ ⁺ + λ ⁻),

{μ} = (1/2) (μ ⁺ + μ ⁻),

where superscripts ± denote left/right values of parameters at an intersected cell's interface.

To illustrate this algorithm we computed two camera views shown in Fig. 1. Computations were performed on 100 × 100 × 200 grid. The length of an edge of each cell is set to 0.1 mm and, therefore, physical dimensions of the computational domain are: (i) 10 mm in x- and y-axis, and (ii) 20 mm in z-axis. The parameter ε was set to 1/2.

In the left Fig. 1a three twisted spirals (the triple helix), which were built from scattering balls, serve as a highly scattering object. The spirals are embedded in a weakly scattering cylinder, whose symmetry axis is aligned along z-axis. It has the background transport coefficient μ = 0.1 mm ⁻¹ and the albedo λ = 0.999. The value of the transport coefficient for each ball is set to μ = 0.75 mm ⁻¹. Two spirals have the background value of the albedo and one absorbing spiral has the value of albedo λ = 0.25. In regions where balls overlap, the density of scattering particles is increased and, therefore, values of parameters are computed accordingly as described in section 2.1. The direct radiation I ₀ enters the domain along the direction s ₀ =2^−1/2 (1,0,−1)^T. The amplitude of the direct radiation is set to 1. A camera rotates around the z-axis. It is convenient to describe a position of the camera by a normal vector n to its CCD array. We set the initial position of the camera as n = (1,0,0)^T . In Fig. 1a the camera was rotated by 153° with respect to the initial position. In the figure on the right (Fig. 1b) two highly scattering cylinders are embedded in a weakly scattering cylinder. The weakly scattering cylinder has the same optical properties as above. Both highly scattering cylinders have μ = 0.75 mm ⁻¹, one of them has a low value of the albedo, λ = 0.25. The direct light enters the domain along the same direction as above. The camera was rotated by 117° from its initial position around the z-axis in the positive direction. In addition to these figures, two multimedia files show animated camera's views over 360° of the triple helix and cylinders (Media 1, Media 2). The camera was rotated around the z-axis with 3° angular step. These two cases will be used below for reconstruction experiments.

2.3. Inverse problem

Depending on the scattering properties of objects the Radon transform [24] can be applied for computing the transport coefficient μ from projection datasets due to domination of the direct radiation in the transmitted light. However, this approach has limited applicability because it does not take into account forward scattered intensity, and does not allow reconstruction of two parameters. Inversion formulae for the attenuated Radon transform ( [25–28]) could not be used with this type of tomographic imaging for the reason that the function B in Eq. (6) depends on the direct radiation, whose direction s0 varies. Moreover, the diffusive nature of light transport in highly scattering regions makes the inverse problem to be three-dimensional. Therefore, we use a variational framework [29].

The variational problem is formulated as a minimization problem of the cost functional:

where the error norm is given as:

 

Fig. 1. (a) (Media 1) Three twisted spirals built from scattering balls are embedded in a weakly scattering cylinder with background transport coefficient μ =0.1 mm ⁻¹ and albedo λ = 0.999. Value of the transport coefficient for each scattering ball is set to μ = 0.75 mm ⁻¹. Two spirals have the background value of the albedo and one absorbing spiral has the value of albedo λ = 0.25. The direct light enters the domain along the direction s ₀ = 2^−1/2 (1,0,−1)^T. Camera was rotated around the weakly scattering cylinder, whose axis is aligned along z-axis, by 153° with respect to the initial position n = (1,0,0)^T, where n is the camera normal. (b) (Media 2) Two highly scattering cylinders are embedded in a weakly scattering cylinder with the same optical properties as in (a). Both highly scattering cylinders have μ = 0.75 mm ⁻¹, one of them has a low value of the albedo, λ = 0.25. The direct light enters the domain along the same direction as in (a).The camera was rotated by 117° from its initial position around z-axis in the positive direction.

Download Full Size | PDF

Intensities I_E and I are experimentally recorded and computed intensities in the direction s, respectively. The function ξ (s) is introduced for convenience and represent sampling of the camera's positions

where −s _n can be replaced with the normal vector to the CCD array n, and N is the number of positions of the camera. Similarly, the functions χ_θ and ς represent sampling of measurements in space and frequency

where M is the number of the camera's pixels; L denotes the number of samples in the Fourier domain (ω); the vector r _m denotes the surface points visible by the CCD camera. Factors σ_m are surface areas supporting rm such that ∫ χ (r)d ³ r gives the total visible area. This form of 𝓔 is chosen in order to simplify a variational procedure. Thus, the function χ allows to replace a sum over surface points visible by the CCD camera with a volume integral. Analogously, the function ς replaces a sum over samples in the Fourier domain with an integral.

The Lagrangian term in Eq. (16) is denoted by 𝓛 and explicitly given by

where 〈·, ·〉 denotes the inner product and J is the adjoint intensity. The dynamic form of the penalty term is chosen, which depends on (k + 1)-th and k-th iterations as

where Δμ = μ _k+1−μ_k, Δλ = λ _k+1−λ_k, α_μ and α_λ are Tikhonov regularization parameters.

The reconstruction algorithm is based on the condition that the first variation δℱ(I,u, J,μ,λ) vanishes. Variation of J recovers Eq. (1) while the variation of I results in the adjoint transport equation

where asterisk denotes complex conjugation. Notice that in Eq. (22) the propagation direction is reversed s = −s _n. Therefore, the adjoint intensity J* propagates from the CCD array in direction of its normal n. Variations of optical parameters μ and λ results in two equations

where functions f_μ and f_λ are computed according to

In many practical cases ω/μc≪1 and, therefore, terms containing this parameter in Eqs. (24) and (25) can be neglected. It is also interesting to notice that the second term in Eq. (24) disappears for the time-independent case. Formally, Eqs. (23) represent an iterative backprojection algorithm, when backprojected functions f_μ and f_λ are products of J* with various combinations of I ₀, I, u and s·∇u.

Tikhonov regularization parameters are estimated from the following consideration. Let us consider α_μ first. From Eqs. (23) we find that

Next, ∥Δμ∥ must decay together with 𝓔^1/2 and, therefore

where the coefficient C ^(k) _μ in general depends on the iteration number k and decreases with k. Analogously, we find

where C ^(k) _λ is analogous to C ^(k) _μ. Iterations are terminated when the functional 𝓔 + ϒ attains its minimal value.

3. Numerical experiments

The dataset for reconstruction was computed by solving direct problems, Eq. (6), when the direction of the incident direct radiation, s ₀, was rotated by 5° around the z-axis over 360° starting from its initial direction s ₀ = (1,0,0)^T. The source plane r ₀ · s ₀ = const was parallel to the z-axis (e _z · s ₀ = 0) and set at some distance from the scattering domain. Free space was assumed to be non-absorbing and non-scattering. For each direction of s ₀ images were acquired by rotating the CCD camera. The camera's viewing direction is controlled by its normal n, which was rotated independently of s ₀ around the z-axis with angular step of 90° starting from the its initial direction n = (1,0,0)^T. The first approximation to the transport coefficient was computed from the projection dataset when s ₀ and n differ by 180°. The albedo was reconstructed from angularly selective intensity measurements, when s ₀ and n differ by 90° and 270°. The first approximation to the transport coefficient was further corrected from angularly selective intensity measurements.

We start reconstruction experiments with the case of two highly scattering cylinders embedded inside a weakly scattering medium. Reconstruction results are shown at z = 10 mm in Fig. 2. The left slice (Fig. 2a) demonstrates reconstruction of the transport coefficient μ, while the right slice (Fig. 2b) displays reconstructed albedo.

For analyzing the three-dimensional structure of reconstructed parameters, isosurfaces were computed and shown in Fig. 3. On the left the isosurface of the transport coefficient is shown (Fig. 3a) and on the right (Fig. 3b) the isosurface of reconstructed albedo is displayed.

As is seen in Fig. 3, the reconstruction is relatively accurate. However, both isosurfaces have kinks along the z-axis. These kinks are caused by obstructions of one cylinder by another for projection acquisition and by the cylinders' shadows for angularly selective measurements. It is worthwhile to mention that obstructions and shadows result in opposite effects on reconstruction.

 

Fig. 2. Reconstruction results showing middle slices at z = 10. (a) Reconstructed transport coefficient μ (b) Reconstructed albedo λ

Download Full Size | PDF

 

Fig. 3. (a) Isosurface of the transport coefficient μ. (b) Isosurface of the albedo λ.

Download Full Size | PDF

 

Fig. 4. Slices showing reconstruction results of the triple helix at two different heights (a–b) Reconstructed transport coefficient μ albedo λ at z = 9mm. (c–d) μ and λ at z = 11mm.

Download Full Size | PDF

Slices showing reconstruction results for the case of a triple helix are displayed in Fig. 4. Slices are taken at different heights. Thus, the first row show parameters at z = 9 mm, while the second row show parameters at z=11 mm. Left column displays the transport coefficient μ and the right column shows the albedo λ. Quantitatively our reconstruction approach demonstrates accurate results. Even increase of the transport coefficient in a regions where scattering balls overlap with each other is clearly seen. However, there is some noise present on the weakly scattering background.

Isosurfaces were also computed for the analysis of spatial accuracy of reconstruction. They are shown in Fig. 5, where (i) Fig. 5a displays the isosurface of the transport coefficient and (ii) Fig. 5b shows the albedo. Reconstruction results demonstrate good separation of the twisted spirals.

Computation time of the direct problem depends on a volume of highly scattering inclusions. That is because the most expensive part of the algorithm is the solution of the Helmholz equation, Eq. (10). For the case presented here the solution time is approximately 29 seconds on 2GHz processor, where 23 seconds are taken for solving the Helmholz equation and only 6 seconds for line integrations, Eqs. (5), (6). The algorithm was implemented in the Fourier domain for further application for time-gating imaging. Thus, for the real case the solution of the Helmholtz equation will be faster. Computation time of the inverse problem depends on the initial guess. Here we assumed a weakly scattering medium initially and, therefore, the first approximation to the transport coefficient was estimated without solving the Helmholz equation. Correction of the transport coefficient and reconstruction of the albedo requires solving the Helmholtz equation. In total, the computation time for solving the inverse problem for each of two cases was less than 2 hours.

 

Fig. 5. (a) Isosurface of the transport coefficient μ. (b) Isosurface of the albedo λ.

Download Full Size | PDF

In summary we have investigated the problem of optical tomography in a weakly scattering medium in the presence of highly scattering inclusions. The approach is based on the assumption that the transport coefficient of the scattering media differ by an order of magnitude for weakly and highly scattering regions. This situation is common for optical imaging of live objects such an embryo and, therefore, we believe that our approach can find applications in biomedical imaging. This approach was verified by reconstructions of two optical parameters from numerically simulated datasets.