Forward problem in diffuse optical imaging

Analytical modelling techniques
Green's functions provide a method for modelling the diffusion equation or the RTE analytically. The Green's function is the solution when the source is a spatial and temporal &delta;-function, from which solutions for extended sources can be derived by convolution. Unfortunately, solutions only exist for simple homogeneous objects (Arridge et al. 1992) or media which include a single spherical perturbation (Boas et al. 1994), although the range of solutions is being extended to more complex geometries such as layered slabs (Martelli et al. 2002).

A number of workers have extended this type of approach. Kim (Kim and Ishimaru 1998, Kim 2004) has devised a method for calculating the Green's function as an expansion into plane-wave solutions, which allows analytical solutions for the diffusion equation or the RTE to be derived. Ripoll et al. (2001) has advocated the use of the Kirchhoff approximation, which models the Green's function between two points in a medium of arbitrary geometry as a sum of the infinite space Green's function plus other Green's functions calculated for diffusive waves which are multiply reflected off the boundary. It can be used alone, or to improve the accuracy of existing modelling techniques. Green's function techniques are now commonly used to solve the forward problem for image reconstruction, particularly for fast imaging techniques where the geometry can be approximated as a slab or an infinite halfspace, such as for optical topography (Culver et al. 2003b, Li et al. 2004).

Spinelli et al. (2003) have developed a perturbation approach first introduced by (Arridge 1995) in which the optical properties are modelled by a Green's function for a slab representing the homogeneous background, with an additional perturbation term representing a spherical insertion. It is well-suited to analysing transmission data measured across a compressed breast with a single, isolated lesion ((Torricelli et al. 2003).

Statistical modelling techniques
Statistical (or stochastic) modelling techniques model individual photon trajectories and have the advantage that the Poisson error is incorporated into the model in a natural and elegant way. The most commonly used statistical technique in diffuse optics, and that which is often regarded as the gold standard to which other techniques are compared, is the Monte Carlo method. The geometry of the model is defined in terms of µa, µs, and the refractive index, and the trajectories of photons, or packets of photons, are followed until they either escape from the object under study or are absorbed. By continuing until the required counting statistics are obtained, data with arbitrarily low statistical errors can be simulated. In optical imaging, Monte Carlo techniques are commonly used to calculate light propagation in non-diffusive domains where the diffusion approximation does not hold (Boas et al. 2002, Okada and Delpy 2003, Hayashi et al. 2003), or to validate results obtained using other, faster, methods (Schweiger et al. 1995, Chernomordik et al. 2002b, Dehghani et al. 2003a).

Random walk theory provides a distinct approach in which photon transport is modelled as a series of steps on a discrete cubic lattice. The time steps may be discrete or continuous (Weiss et al. 1998). Random walk theory is particularly suited to modelling time-domain measurements (Chernomordik et al. 2000) and has been used, for example, to determine the time spent by photons in a scattering inclusion (Chernomordik et al. 2002a) and to quantify the optical properties of a breast tumour (Chernomordik et al. 2002b). Recently, an extension to the technique has been developed for modelling media with anisotropic optical properties (see section 3.3.6), maintaining the cubic lattice but allowing the transition properties along different axes to differ (Dagdug et al. 2003). Similarly, Carminati et al. (2004) have developed a more general method based on random walk theory which allows the transition from single-scattering to diffusive regimes to be explored.

Numerical modelling techniques
Numerical techniques are required if more complex geometries are to be modelled. The natural choice for representing the inhomogeneous distribution of optical properties in an arbitrary geometry is the finite element method (FEM) which was first introduced into optical tomography by Arridge et al. (1993). The method is explained in detail by Arridge and Schweiger (1995), Arridge (1999) and Arridge et al. (2000). FEM has become the method of choice for modelling complex inhomogeneous domains in optical imaging (Bluestone et al. 2001, Dehghani et al. 2003b), although the finite difference method (FDM) (Culver et al. 2003a, Hielscher et al. 2004), finite volume method (FVM) (Ren et al. 2004) and boundary element method (BEM) (Ripoll and Ntziachristos 2003, Heino et al. 2003) have been used in more specialised applications.

The finite element method requires that the reconstruction domain be divided into a finite element mesh. In principle, this is a completely general technique which can be applied to any geometry. In practise, however, it is difficult to create a finite element mesh of irregular objects with complex internal structure, and the development of robust, efficient 3D meshing techniques is still the subject of active research. A particular challenge is meshing the head, while respecting the convoluted internal boundaries of the scalp, skull, cerebrospinal fluid (CSF), grey matter and white matter derived from a segmented MRI image. Inclusion of such anatomical details has been shown to improve the quality of reconstructed images in EIT (Bagshaw et al. 2003), although generating the 3D finite element mesh was difficult and time-consuming (Bayford et al. 2001, Lionheart 2004). The use of realistic finite element meshes in optical imaging has progressed from 2D models of internal structure taken from segmented MR images of the head and breast (Schweiger and Arridge 1999, Brooksby et al. 2003), through 3D finite element meshes with complex surface shape but no internal structure (Bluestone et al. 2001, Gibson et al. 2003, Dehghani et al. 2004), to a recent report by (Schweiger et al. 2003) of simulations from 3D finite element meshes of the breast and head with anatomically realistic internal structure. A further complication is the consideration that optimal computational efficiency requires a finite element mesh which can represent the internal field adequately whilst using the smallest possible number of elements. One approach is to adaptively refine the mesh, placing more elements where the field changes most rapidly (Molinari et al. 2002, Joshi et al. 2004). Another approach, which is suited to modelling a segmented volume taken from MRI, is to resample the pixels in the MRI image onto a regular grid which can then be solved with FDM (Barnett et al. 2003)

Use of prior information
Perhaps the most significant recent trend in optical imaging has been the inclusion of prior anatomical information into the forward problem, most commonly in the form of an anatomical MRI image. Prior information was first used in optical tomography reconstruction by Pogue and Paulsen (1998), and by Schweiger and Arridge (1999). The latter used a more sophisticated reconstruction procedure, which began by building a finite element mesh with four regions segmented from an MR image. They smoothed the model to account for gradual transitions between regions and generated simulated data, to which noise was added. They then used a two-step reconstruction process in which the optical properties of the regions were reconstructed first, followed by a full reconstruction onto the finite element mesh basis. The use of prior anatomical information was shown to improve the image quality in all cases. Ideally, the anatomical prior and the optical image should be acquired simultaneously so the two images are correctly registered. This impacts on the design of the image acquisition system and the experimental procedure and so these approaches are reviewed in the appropriate section on clinical applications below.

Non-diffusive regions
Light transport can be modelled adequately using the diffusion approximation in most biological tissues. Among the exceptions to this are tissue volumes which are smaller than a few scattering lengths, and regions where µa is comparable to or greater than µ's, such as the CSF which surrounds the brain and fills the central ventricles. These regimes are generally modelled using higher order approximations to the RTE.

It is generally too computationally expensive to solve the RTE fully in a practical reconstruction scheme. Instead, approximations are made such as the method of discrete ordinates, which solves the full RTE on a regular grid using FDM or FVM by quantising the allowed directions of travel, $$\hat s$$. This makes the problem tractable but because only preferred angles of travel are allowed, there may be regions where light cannot reach. This is known as the Ray Effect and restricts the application of the method. The method of discrete ordinates was applied to optical imaging by Dorn (1998) and has been used successfully (Klose and Hielscher 1999, Ren et al. 2004) for imaging small objects such as the finger (Hielscher et al. 2004), and in fluorescence and molecular imaging of small animals (Klose et al. 2004). Alternatively, the RTE can be solved using the PN approximation (Aydin et al. 2002) or by expansion into a rotated spherical harmonic basis (Markel 2004).

A different approach can be taken for modelling light transport in the head, where the healthy CSF is non-scattering (although brain injury may cause proteins and blood products to leak into the CSF, making it diffusive (Seehusen et al. 2003)). If an anatomical MRI is available, then it is possible to segment the head into diffusive and non-diffusive regions. A coupled radiosity-diffusion model has been developed at UCL which models light transport in scattering regions using the diffusion model, and in clear regions using a visibility model which couples points on the void boundary which are mutually visible. The model has been applied to 2D (see Figure 2) and some 3D geometries (Riley et al. 2000), but it has not yet been applied to a 3D head-like model. A similar approach has been adopted by Schulz et al. (2003, 2004) who use a coupled model for non-contact molecular imaging where the animal is modelled diffusively and the space between the optical fibres and the animal is modelled using a free-space propagation model. One disadvantage of the current implementation of the radiosity-diffusion model is that the boundary of the clear region must be known.

Anisotropy
Another regime in which the diffusion approximation does not apply is that of anisotropic scattering. Certain tissues such as the skin, nervous tissue and muscle cells are anatomically anisotropic at the cellular level, and this leads to anisotropic scattering properties (Nickell et al. 2000).

Heino and Somersalo (2002) derived a modified form of the diffusion equation in which the diffusion coefficient, $$\kappa=1/{3(\mu_a+\mu_s[1-g])}$$, is replaced by a diffusion tensor $$\kappa=1/{3([1\mu_a+\mu_s]1-\mu_sS)}$$, where S is (in 3D) a 3×3 matrix whose diagonal contains the anisotropic scattering coefficients. They make the anatomically realistic assumption that µa is isotropic everywhere and distinguish between knowing a priori the direction but not the strength of the scatter anisotropy (as may be obtained from MR diffusion tensor imaging (Le Bihan et al. 2001)), and having full knowledge of both the direction and the strength of the anisotropy. Such prior information is required if the solution is to be unique. Statistical reconstructions of simulated data demonstrated that artefacts in µa will result if anisotropic µ's is not fully considered. They validated the results by comparison with a Monte Carlo model (Heino et al. 2003).

Dagdug et al. (2003) studied the same problem using a random walk formulation in which the transition probability is allowed to vary with direction. To date, this method only allows the direction of anisotropy to be parallel to the co-ordinate axes. The model has been shown agree with time-resolved measurements on anisotropic phantoms (Hebden et al. 2004).