1 Introduction

Optical Tomography is an example of a non-linear illposed inverse problem. A complete treatment of image reconstruction therefore requires either implicitly or explicitly the optimization of an objective function with respect to the parameters of a model[1]. The most commonly used model is the diffusion approximation which is the lowest order non-trivial approximation to the more correct Radiative Transfer Equation (RTE)[1, 2]. The reason for its success is that tissues of interest in clinical applications, such as the human breast and peripheral muscles, are of a scale such that detected photons have undergone a high enough number of scattering events that the dependence of fluence on direction is only linear.

However, one major application interest for optical tomography is its use for imaging of the brain. Within the head the regions of cerebrospinal fluid (CSF) around the brain and in the ventricles are non-scattering while still absorbing. The presence of CSF prevents the accurate modeling of photon propagation within the regions of interest when using the diffusion approximation and leads naturally to the attempt to develop a modelling scheme based on the RTE[3, 4, 5, 6]. The principle difficulty then is that accurate handling of the non-scattering regions implies the use of a very computationally expensive numerical method.

One proposed alternative to the full RTE is the Radiosity-Diffusion model which assumes diffusive regions coupled by non-scattering voids [7, 8]. Propagation in the void is handled via geometrical optics and the interface requires the development of non-local boundary conditions [9]. The advantage is the greatly increased computational efficiency whilst maintaining good accuracy. This efficient model has allowed the investigation of the effect of voids on the accuracy of image reconstruction in 2D [10]. However optical tomography is in reality a 3D problem [11]. In this paper we report the development of the radiosity-diffusion model in 3D and present the first assessment of optical tomography in a 3D geometry involving voids.

2 The Radiosity-Diffusion Model

Let Ω represent a domain consisting of R diffusing regions {Ω₁, Ω₂, … Ω _R } and V void regions {Ξ₁, Ξ₂,… Ξ _V }. Let Ω ^d =${\cup }_{i=1}^R$Ω _i be the union of all diffusing regions and Ξ ^d =${\cup }_{i=1}^V$ Ξ _i be the union of all void regions. Thus Ω=Ω ^d ∪ Ξ ^d . Each region has an outer boundary ${\mathit{\partial }}_i^{+}$ and a non-negative number of inner boundaries ${\mathit{\partial }}_{i\mathit{,}k}^{-}$. We define ∂Ω _i =∂${\mathrm{\Omega }}_i^{+}$ ∪(∪ _k ∂${\mathrm{\Omega }}_{i\mathit{,}k}^{-}$ ) for diffusing regions, with similarly ∂Ξ _i =∂${\mathrm{\Xi }}_i^{+}$ ∪(∪ _k ∂${\mathrm{\Xi }}_{i\mathit{,}k}^{-}$ ) for voids. The outer boundary ∂${\mathrm{\Omega }}_1^{+}$ will also be denoted by ∂Ω.

We consider a frequency-domain system so that the photon density Φ(r; ω) at modulation frequency ω satisfies the homogeneous equation

where κ=1/(3(µ _a+µ′ _s)) is the diffusion coefficient defined in terms of µ _a(r) and µ′ _s(r), the spatially varying absorption and reduced scattering coefficients respectively and with local inhomogeneous boundary conditions

where η is the incoming flux modelled as a Neumann source [13], and v is the surface normal pointing into the void region, and non-local boundary conditions

and h(m, m′) is unity if m, m′ are in line of sight across the void, and zero otherwise. A more exact boundary condition replaces the term Φ(m′;ω)/2A in eq(3)by

$\left(\Phi (m^{\prime };\omega )-2\kappa \left(m^{\prime }\right)\frac{1+R^{\left(1\right)}}{1-R^{\left(0\right)}}\frac{\partial \Phi (m^{\prime };\omega )}{\partial \nu }\right)(1-{\mid R\left(\theta \right)\mid }^2)$

where R ⁽⁰⁾, R ⁽¹⁾ are derived from the Fresnel coefficients taken over local coordinates

R ⁽⁰⁾=2${\int }_{\mathrm{\pi }/2}^{\mathrm{\pi }}$ R(ϑ)sinϑcosϑdϑ, R ⁽¹⁾=3${\int }_{\mathrm{\pi }/2}^{\mathrm{\pi }}$ R(ϑ)sin ϑcos²ϑdϑ

See [12] for a detailed derivation. These two forms were compared in [9] and found to give comparable results. The constant for the domain boundary in eq(2) and (3) is given by A=(1-R ⁽¹⁾)/(1-R ⁽⁰⁾). For each Neumann source term η the measureable is

3 Inverse Problem

In the adjoint field approach [14] we assume a finite number of input fluxes {η_j ; j=1,…, S}, and given data {g_j ;j=1,…, S}, and look for the minimisation of the norm

A descent direction for minimising C with respect to the data from the j^th source is given by

with ψ _j the solution to the adjoint diffusion equation

with boundary conditions

For the optimization scheme we employ the same conjugate-gradient scheme as used in [10]. However, adequate treatment of the 3D problem requires some special attention, which we describe in the next section.

4 Implementation in 3D

Our numerical method utilises the Finite Element Method (FEM) to allow application to general domain boundaries and parameter distributions. In the results presented here we use linear tetrahedral elements generated using the public domain package Netgen [15], which implies that the fluence Φ is expressed as a linear combination of piecewise linear shape functions {u_k ; k=1… D} with value one at nodal points N_k and zero at all other nodes. For simplicity we consider the case of concentric spherical surfaces which are defined parametrically. The discretisation of eq(1) leads to sparse symmetric matrices, whilst the implementation of the non-local boundary conditions eq(3) leads to a dense coupling matrix whose order, D_B , is the number of nodes in the void boundary. In 3D this number is proportional to the surface area of the void surfaces which is the principle extra overhead to the computational requirements of the foward solver. In the previous 2D treatment [8] we calculated the visibility function by a full search of all void nodes which becomes computationally expensive in 3D. For the case considered here we can determine the visibility function h(m, m′) as well as the surface normal functions $\widehat{\nu }$ directly from the parametric representation of the surfaces which greatly simplifies the setup time of the forward solver. The function

extracted from eq(3), is called the form factor corresponding to its use in Computer Graphics [16]. In our current implementation we approximate this function by a bilinear expansion over the FEM shape functions on the surface of the void [17]. Between two element faces τ_α , τ_α′ ; ∊ δΞ containing N, N′ nodes respectively the form factor can be approximated by

where n(k) maps the local node to the global nodes in the complete mesh.

5 Results

In fig.1 we show the geometry being considered. We take a sphere radius 25mm with concentric sphere radius r_inner and a concentric void gap with outer radius r_outer . 32 sources and detectors are arranged in three rings. We used background parameters µ _a=0.01 mm ^-1, µ′ _s=1 mm ^-1 with µ _a=0.005 in the void region. For the sources and detectors we use a cosine weighted patch in the parameters of the surface representation. For the reconstruction basis we use (20×20×20) tricubic interpolated voxels.

 

Fig. 1. Left : cutaway of spherical mesh. Right : location of sources and detectors on the sphere surface

Download Full Size | PDF

 

Fig. 2. distribution of photon density, photons/mm³ (top row) and mean time of photon flight, picoseconds (bottom row) over sphere surface. Left to right, solid sphere (no gap), gap widths 3mm, 4mm, 5mm.

Download Full Size | PDF

In fig.2 we show the surface data generated by a single source for the case of a solid sphere and for three gaps defined by r_inner =17,16,15mm, and in each case r_outer =20mm. The difference between each void case and the solid sphere emphasises the increased light intensity and decreased mean time in the gap region.

 

Fig. 3. Sensitivity functions for the 3mm gap case. Left intensity (photons/mm²), right mean time(picoseconds mm). The functions plotted are cross-sections through the equatorial plane of the sphere. Also available as a QuickTime movie, pmdf.mov. (3.8MB)

Download Full Size | PDF

In the animation attached to fig.3 we show the sensitivity functions in absorption for intensity (y(ω=0)) and mean time $\frac{\partial \ln \mathit{y}}{\partial \omega }{\mid }_{\omega =0}$ measurement types, in the case of a 3 mm wide gap, for 16 different detector locations around the equator of the sphere. The sensitivity functions are derived by using eq(6) with ψ _j an adjoint Green’s function obtained using a δ-function for the local boundary condition in eq(8); see [1, 18] for a more detailed derivation, and the related paper [6]. For intensity measurements the sensitivity is mainly confined to the regions around the source and detector sites in the outer ring, whereas for mean time the gap has less of an adverse effect, providing higher sensitivity in the inner region. We use mean time data only in the following for consistency with [10].

 

Fig. 4. Target images (top row) and reconstructions (bottom row) for the 3mm gap case. The images are transverse, sagittal and coronal slices through the true centre of the blob, orientated according to the diagram in the top right panel. Bottom right shows a profile along the equatorial diameter through the blob centre. A movie showing a rotating orthographic view is attached rotating3d.mov (3.8MB).

Download Full Size | PDF

In fig.4 we show the reconstruction of a blob in the 3mm gap. The blob had a radius of 3mm and an absorption of 0.02mm ^-1 with its center placed at position (12,0,0). The images shown are at iteration 40 of the conjugate gradient scheme. The reconstruction shows good localisation despite quite a large degree of broadening.

6 Discussion and conclusions

The correct treatment of a non-scattering region in diffusing media is an important requirement for application of Optical Tomography to the brain. In this paper we showed initial results of a 3D extension to the Radiosity-Diffusion model that provides one approach to this problem. Even though the cases considered here are simple, they indicate that, if the location of the void is known, the signal detected on the domain surface still carries enough information to allow a reconstruction. A more complete treatment should consider data types in addition to just the mean time as used here, as well as simultaneous reconstruction of absorption and scattering. Furthermore, one of the interesting results found in [10] was the improvement in reconstruction quality when the boundary of the void gap is not smooth. Investigation of this case requires a more sophisticated determination of form factors which is under development.

In the case where the void location is not known, we are investigating the use of a boundary recovery scheme related to the method reported in [19]. The application of this scheme to the Radiosity-Diffusion model requires the determination of differentiability of the forward operator at the diffuse/non-diffuse interface.