Optical imaging/recon

Background
The 1997 review of modelling and image reconstruction techniques (Arridge and Hebden 1997) is still a useful introduction to the basic theoretical concepts. Progress since 1997 has largely focussed on developing more realistic and efficient models of light transport in tissue, and on solving the ill-posed inverse problem in an increasingly rigorous way. In particular, it is becoming increasingly common to include prior information of the anatomy and optical properties of tissue in both the modelling and the image reconstruction.

Three processes are required to reconstruct an optical image. First, the transport of light in tissue must be modelled. Second, this model must be used to predict the distribution of light in the object under examination. This stage, the forward problem, allows the measurements to be simulated from the model, and generates (explicitly or implicitly) a sensitivity matrix (the Jacobian of the forward mapping) which relates the measurements to the internal optical properties. Finally, the image is reconstructed by inverting the Jacobian and solving the inverse problem. In x-ray CT, scatter is minimal, so the forward problem becomes a series of integrals along lines connecting the sources and detectors (a Radon transform), and the inverse problem is linear and well-posed. In optical imaging, however, scatter is dominant, meaning that the forward problem becomes a series of integrals over the entire volume. Each measurement is therefore sensitive to the whole volume, meaning that the inverse problem is ill-posed and underdetermined and that complex, computer-intensive reconstruction techniques are required. For recent reviews of image reconstruction in optical imaging, see Arridge (1999), Kohlemainen (2001) and Schweiger et al. (2003).

It should be noted that similar issues are being addressed in electrical impedance tomography (EIT), which is a related medical imaging technique in which surface measurements of electrical impedance are used to reconstruct volume images of tissue conductivity (Metherall et al. 1996, Borcea 2002, and Lionheart 2004).

Modelling light transport in tissue
Before image reconstruction can be attempted, a model of photon transport in tissue is required. For diffuse optical imaging, the diffusion equation, in which the propagation of light is assumed to be isotropic, is commonly used. In this section, we discuss the derivation of the diffusion equation from the radiative transport equation. Arridge (1999) provides a more detailed review.

A full description of light propagation in tissue is provided by the radiative transport equation (RTE; Equation 1). This is a conservation equation which states that the radiance (the number of photons per unit volume),, for photons travelling from point r in direction  at time t is equal to the sum of all the mechanisms which increase  minus those effects which reduce it. Equation (1) shows the time-domain RTE. The RTE in the frequency-domain is obtained by replacing  by.

(1)

In (1),  is the scatter phase function, which gives the probability of a photon scattering from direction  to,   and   is the light source at r at time t travelling in direction. For clarity, we use the same symbols as Arridge (1999).

The RTE is an approximation to Maxwell’s equations and has been used successfully to model light transport in diffusive media, turbulence in the earth’s atmosphere, and neutron transport. However, it does not include wave effects, so the wavelength must be much smaller than the dimensions of the object under study. It also requires the refractive index to be constant in the medium, although extensions for spatially varying refractive index have been derived (Ferwerda 1999, Marti-Lopez et al. 2003, Tualle and Tinet 2003). Unfortunately, to solve for the properties of the light field at all points in a large 3D volume, as is required in optical tomography, the RTE is extremely computationally expensive and simpler models are generally sought. Kim and Keller (2003) recommend the use of a modified Fokker-Planck equation as an approximation to the RTE where scatter is strongly forward-biased.

Three variables in the RTE depend on direction : the specific intensity ?, the phase function, and the source term q. If these are expanded into spherical harmonics, we obtain an infinite series of equations which approximate to the RTE. The PN approximation is obtained by taking the first N spherical harmonics, which gives (N+1)2 coupled partial differential equations. As N increases, the PN approximation models the RTE more accurately, but with increasing computational requirements (Aydin et al. 2002).

If we now take the P1 approximation and assume that the phase function  is independent of the absolute angle, i.e. , that the photon flux changes slowly and that all sources are isotropic, we obtain the diffusion equation (2):

,		(2)

where photon density, and the diffusion coefficient .

The diffusion equation (2) has been widely and successfully used to model light transport in tissue, although it is necessary to assume that light propagates diffusively. This is generally the case in bulk tissue but the assumption breaks down near the source, near the surface, near internal boundaries, in anisotropic tissues, and in regions of either high absorption or low scatter. In these situations, higher order approximations to the RTE may be required.

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 ?-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. 2003c, 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. (2000a). 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 1999b, 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 (1999b). 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,. 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,, is replaced by a diffusion tensor, 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. 2004a).

Linear image reconstruction
The forward problem, discussed above, involves calculating simulated data y (which may be CW, time-domain or frequency-domain data), given a forward operator F and knowledge of the sources q and internal optical properties x (which may include µa and µ's). The forward problem may therefore be formulated as:

.						(3)

To reconstruct an image, it is necessary to solve the inverse problem, which is to calculate the internal optical properties x, given data y and sources q (Arridge 1999):

.						(4)

This is a non-linear problem but it can be linearised if the actual optical properties x are close to an initial estimate x0 and the measured data y is close to the simulated measurements y0. This is typically the case in difference imaging where measurements are taken before and after a small change in the optical properties. Then we can expand (3) about x0 in a Taylor series:

(5)

where and  are the first and second-order Fréchet derivatives of F, respectively. The Fréchet derivative is a linear integral operator mapping functions in the image space to functions in the data space. In some cases, such as in section 3.3.1, the kernel of the integral operator is known in terms of the Green’s functions. This approach was taken by Boas et al. (1994).

More generally, the forward problem can be solved by a numeric technique by representing the Fréchet derivatives and   by matrices J (the Jacobian) and H (the Hessian), respectively. Equation 5 can then be linearised by neglecting higher order terms and considering changes in the optical properties  and data to give the linear problem (6):

.						(6)

Linearising the change in intensity in this way gives rise to the Born approximation; linearising the change in log intensity instead is the Rytov approximation and may lead to improved images by reducing the dynamic range of y. Either way, image reconstruction consists of the problem of inverting the matrix J, or some form of normalised J. This may be large, underdetermined and ill-posed but standard matrix inversion methods can be used. These techniques differ in the way in which the matrix inversion is regularised to suppress the effect of measurement noise and modelling errors. Perhaps the most common techniques are truncated singular value decomposition, Tikhonov regularisation, and the algebraic reconstruction technique (ART) (Gaudette et al. 2000). The Moore-Penrose inverse  offers a more efficient inversion if J is underdetermined and leads naturally to a Tikhonov- type formulation (7) where I is the identity matrix and ? is a regularisation parameter. J itself is calculated from the forward model, most efficiently using the adjoint method of Arridge and Schweiger (1995). The linear inverse problem can then be expressed as:

.				(7)

Non-linear image reconstruction
If the inverse problem cannot be assumed to be that of reconstructing the difference between two similar states, but is instead a reconstruction from a single acquisition to obtain absolute values of the optical properties, the linear approximation cannot be justified and the full non-linear problem must be solved. This is commonly the case in breast imaging (Dehghani et al. 2003b), and in static imaging of the brain (Bluestone et al. 2001, Hintz et al. 2001, Hebden et al. 2002). To solve the non-linear problem, we define an objective function, which represents the difference between the measured data y and data which is simulated using the forward model F(x). If  is the distribution of optical parameters which minimises, then it is also the model which best fits the data, and is therefore taken to be the image. However, because the problem is ill-posed, it must be regularised, and can be expressed as:

,				(8)

where  is a regularisation parameter, and   represents prior information. This may be very simple, for example, in which case (8) becomes an iterative solution of (7). However, it is increasingly common for  to include anatomical and other information, see section 3.4.3. represents the L2-norm and gives the least squares solution. Equation 8 can be extended to include the covariance of the data CP and the covariance of the image CQ as follows:

.				(9)

This approach gives a weighted least squares solution which reduces the influence of noise and cross-talk both in the data and in the image. Measuring CP and CQ may not be straightforward. Equation 9 is a non-linear minimisation problem which is generally solved either by a Newton method such as the Levenberg-Marqardt algorithm, or by a gradient method such as conjugate gradients (Arridge 1999).

Uniqueness
No general uniqueness results exist for optical imaging, because of the range of different unknowns (µa, µ's, refractive index, coupling coefficients, anisotropy, geometry, etc.) and the range of different measurables (time-domain or frequency- domain measurements or CW measurements of intensity alone). Some useful results do exist, however, to illustrate situations which may not yield unique solutions (Isakov 1993).

The first important uniqueness result directly relevant to the diffusion equation was derived by Arridge and Lionheart (1998) and developed further by Arridge (1999). They demonstrated that measurements of intensity alone cannot separate the effects of absorption and scatter. Furthermore, if the refractive index is also unknown, even frequency-domain or time-domain measurements are not sufficient for uniqueness, although Matcher (1999) showed that uniqueness may be restored if the underlying model is taken to be the P1 approximation. The specific case of optical mammography using a slab geometry was analysed by Romanov and He (2000), who showed that if full time-domain data is available then reflection measurements (where the sources and detectors are on the same side of the slab) can uniquely determine a piecewise continuous distribution of either µa or µ's if the other parameter is given, but transmission measurements can distinguish µa and µ's if both are unknown. In an interesting application of uniqueness theory, Corlu et al. (2003) used a uniqueness argument to determine which wavelengths gave the optimal separation between different chromophores in spectroscopic imaging, and to reconstruct directly for chromophore concentration.

Arridge and Lionheart (1998) illustrated their uniqueness result with an example of two distributions of µa and µ's which gave indistinguishable intensity data on the boundary, and Hoenders (1997) described a whole class of solutions which could not be distinguished from one another. More specifically, Schweiger and Arridge (1999a) plotted the objective function of the inverse problem against µa and µ's for a range of different datatypes extracted from the TPSF and showed that, where certain combinations of datatypes gave a well-behaved minimum, measurements of intensity alone gave an objective function with an extended minimum at approximately µa µ's = constant, suggesting that µa and µ's could not be separated.

However, as Hoenders (1997) points out, the existence of a solution to an inverse problem depends on more than just the uniqueness. It also depends on other issues including the quality of the data and the statistical behaviour of the problem. In particular, even if a problem is non-unique, it may be possible to distinguish between two identical data sets if appropriate prior information is available, whether it is used implicitly in the form of regularisation or normalisation of the Jacobian (Pei et al. 2001, Xu et al. 2002a), or explicitly. This may enable the theoretical results outlined above to be reconciled with experimental work by (Schmitz et al. 2002) and (Pei et al. 2001), and by (Jiang et al. 2002) and (Xu et al. 2002a) who have reconstructed images showing some separation between µa and µ's from intensity data alone. It remains to be seen how non-uniqueness affects reconstruction of clinical images, where the potential for a reconstruction to converge to an incorrect solution may be unacceptable.

Use of prior information
Advances in modelling and reconstruction techniques cannot overcome the fact that diffuse optical imaging is a non-unique, ill-posed, underdetermined problem and that this puts a limit on the image quality, particularly the spatial resolution, which can be achieved. One way to improve the quality of the image reconstruction is to make maximum use of prior information, which can be obtained from anatomical imaging techniques or by considering the physics and the physiology of the problem. This reduces the effect of the ill-posedness by improving the accuracy of the model, and makes the problem better determined by allowing the small number of measurements which are obtained to be used in a more effective way. For example, Ntziachristos et al. (2000) used an MR image of a breast to provide the location of a tumour a priori. Using this approach, the spatial resolution of the optical image effectively becomes that of the MR image. The reduction in quantitative accuracy due to the partial volume effect seen in optical images reconstructed without prior information does not occur, and truly quantitative images of the haemodynamic properties of the tumour can be obtained.

The most straightforward way to include prior information into the image reconstruction is to use an anatomically realistic forward model (section 3.3.3). This increases the accuracy of the forward model so that it can represent the measurements more precisely, thereby improving the image reconstruction. Further improvements can be made by including information about the covariance of the data CP, and of the image CQ, as in (9). The diagonal of CP is simply the variance of each measurement, and gives an indication of the reliability of that measurement. Off-diagonal elements are the covariance between each pair of measurements and give an indication of the interdependence between measurements. Introducing CP into (9) transforms the fit between the data and the model into a basis where all measurements are independent. CQ contains information about the predicted smoothness of the image. Prior knowledge of the anatomy can also be incorporated into the inverse problem. In the example of a breast tumour given above, the change in optical properties is assumed to come from either a region of interest defined from anatomy (Ntziachristos et al. 2002c) or be heavily biased towards that region (Brooksby et al. 2003, Li et al. 2003). The degradation in image quality due to uncertainties in the anatomical prior, and the best way to minimise this degradation, are as yet uncertain (Schweiger and Arridge 1999b). Including prior anatomical information has advantages for both linear and non-linear reconstruction (Boas et al. 2004a)

Another approach is to express (8) as a statistical inverse problem in a Bayesian framework. Bayes’ theorem can be written as  where is the posterior probability of obtaining the image x given the data y,  is the likelihood of obtaining y given x (which is obtained by solving the forward problem), and  is the prior density, which in these terms is simply an estimate, made before carrying out the experiment, of the most likely image. If we assume that the distribution of errors is Gaussian, then Bayes’ theorem takes the form

(10)

where the likelihood  and the prior. The value of x which maximises (10) is the maximum a posteriori (MAP) estimate. Maximising (10) is equivalent to minimising its negative logarithm, which is (8), but can offer further approaches which deal with prior information and covariance estimates in a direct and intuitive manner (Kwee 1999, Kohlemainen 2001, Mosegaard and Sambridge 2002, Evans and Stark 2002, Oh et al. 2002).

Recently, efforts have been made to incorporate Monte Carlo techniques into the inverse problem. Standard Monte Carlo techniques are too slow to be used for practical image reconstruction, although faster methods are available. One example is the Metropolis-Hastings algorithm which is an example of a Markov Chain Monte Carlo (MCMC) technique, in which individual photon steps are drawn from a biased random walk which is chosen to minimise the number of discarded photons (Mosegaard and Sambridge 2002). These techniques have been used as the basis of image reconstruction algorithms in EIT (Kaipio et al. 2000, West et al. 2004), partly because the Monte Carlo simulations form a probability density which can be interpreted as the Bayesian maximum a posteriori estimate, leading very naturally to a statistical expression of the inverse problem (Evans and Stark 2002).

The approaches discussed above assume that the prior information, whether in terms of the data covariance or the image statistics, is governed by Gaussian statistics. This is not necessarily the case: a Poisson model for the data may be more appropriate if photon counting errors are significant, and assuming that the pixel values are distributed as Gaussian random variables forces the image to be homogeneously and isotropically smooth. The latter issue has been identified as a significant weakness and has been addressed in two ways. First, Kaipio et al. (1999) introduced an inhomogeneous, anisotropic image covariance term which allowed, for example, a change in optical properties in the brain to be correlated more closely with neighbouring pixels in the brain than with those in the skull. A second approach is to minimise the total variation, defined as the integral of the absolute gradient of the image, rather than the L2-norm (Borsic 2002). The L2-norm forces the image to be smooth, and is optimal if the optical properties are distributed as a Gaussian random field. If the distribution of optical properties is known to be piecewise constant, the L2-norm will smooth the image, whereas minimising total variation will reduce oscillations in the image but still allow boundaries to be sharp. Other non-Gaussian priors may be considered.

Methods that exploit symmetry
Some improvements to speed and robustness of the inverse problem can be obtained for systems that possess symmetry. Examples are rotational symmetry in a cylinder or sphere, or translational symmetry in a slab. It is easiest to develop this idea for the analytical version of the linear inversion kernels based on Green's functions. Green's functions are the kernels of the inverse solver for a PDE, and when the domain possesses symmetry, the PDE may often be solved by a separation of variables technique, which implies that the inversion kernels can be expressed as the product of functions in the separated variables. In the work of Markel and Schotland (2001), the slab geometry gives rise to a product of transverse plane waves and one-dimensional integral equations in the perpendicular direction. It is this decoupling of transverse propagating waves, and perpendicular scalar waves that allows faster inversion of a set of one dimensional integral equations, rather than one coupled 3D integral equation. The solution of these linear equations builds up a solution in the basis of the transverse component, which can then be transformed into image space by a fast transform method. Similar methods apply in the cylindrical and spherical geometries where the transverse components are Fourier and spherical harmonic waves respectively. Furthermore the case where only part of the data is sampled can be included by convolution with a masking function. Finally Markel et al. (2003) considered extension to the nonlinear case by applying an iterative method, linearised around the symmetric case.

As well as the analytical case, a discrete implementation such as FEM can also be developed. This approach was developed by Metherall et al. (1996) for EIT, and applied to optical tomography by Hampel and Freyer (1998).

Solving for coupling coefficients
The measured data in optical tomography depends not only on the properties of the object under examination, but also on the characteristics of the source and detector fibres, instabilities in the source power and detector efficiency, and on the efficiency with which light is coupled into and out of the medium. A well-designed calibration procedure can minimise the effects of the optical fibres and the instrumentation (Hillman et al. 2000, McBride et al. 2001) but characterising the coupling between the fibres and the tissue can be more difficult. One approach which has been successfully adopted, particularly with regard to optical topography, is difference imaging, in which images are reconstructed using the difference between measurements recorded at two states in such a way that coupling effects cancel (Section 4.1). Alternatively, images can be reconstructed using temporal data only (Hebden et al. 2002), although this may reduce the ability to separate the absorption and scattering properties. These approaches are not always appropriate, however, and so techniques have been developed to correct for coupling effects during image reconstruction.

Schmitz et al. (2000) and Boas et al. (2001b) introduced a method for solving for the unknown coupling coefficients as part of the image reconstruction problem. Each source and detector was associated with a complex coupling term. A measurement Mij between source i and detector j was modelled as  where Si and Dj are the source and detector scaling factors and mij is the “ideal” measurement, depending only on the internal optical properties of the object. The reconstruction problem is then posed in terms of minimising the difference between the measured data and , where Si and Dj are solved for in the inverse problem. Images could be reconstructed successfully if up to 80% uncertainty in amplitude was added to simulated data (Boas et al. 2001b). The technique has since been extended to allow small errors in optode position to be corrected (Stott et al. 2003, Culver et al. 2003c) and has been used for dynamic breast imaging (Intes et al. 2003). Vilhunen et al. (2004) have implemented a related method which can be used when the sources and detectors have rotational symmetry. Oh et al. (2002) used a Bayesian formulation to reconstruct 3D images non-linearly by treating the unknown source-detector coupling coefficients as unwanted “nuisance” parameters. 3.4.7	Use of dynamic information Image reconstruction is generally applied to the spatial domain but additional information can be obtained by examining the temporal evolution of the signal. Physiological signals from the heart rate, ventilation and blood vessels are time- dependent. If the image acquisition time is short compared to the physiological fluctuations, an image time series can be obtained, which can then be processed to give dynamic images showing parameters such as the covariance or the amplitude and phase at a particular frequency of physiological interest (Barbour et al. 2001, Graber et al. 2002).

A different situation occurs when the image acquisition rate is not short compared to the physiological changes. In this case, the optical properties change during the acquisition of each image and so a straightforward static image reconstruction would not be valid. The Kalman filter technique has been used to model the image space as a state whose properties evolve with time in a known manner, the details of which are supplied as prior information. It may be expressed as:

,				(11)

where xt is the image x at time t, K is the state transition matrix and nt is additive noise at time t. In the simplest case, K is the identity matrix, in which case the update is a random walk. This approach helps to condition the inverse problem by using the current state, and knowledge of how it evolves, to constrain the possible solutions for the next state. It can be used to image events which occur on a shorter timescale than the image acquisition by calculating an updated image using the Kalman filter when only a limited subset of data has been acquired. This is ideally suited for systems which acquire data from sources sequentially but the detectors simultaneously, in which case an updated state can be estimated when each source is activated. This approach has been successfully applied to both simulated data (Kohlemainen et al. 2003) and real data (Prince et al. 2003) in optical tomography. Eppstein et al. (2001, 2002) have developed an image reconstruction technique for fluorescence tomography which uses the Kalman filter to update each iteration of a static non- linear reconstruction procedure.

Recovery of object shape
Another approach to optimising the use of data in optical tomography is not to reconstruct for a large number of image pixels, which leads to a highly underdetermined problem, but instead to reconstruct for the boundaries and properties of internal regions which can be assumed to have piecewise constant optical properties.

One method to solve for the boundaries of internal regions is to formulate a minimisation problem where the image parameters are the Fourier coefficients of the smooth region boundaries. This was first solved assuming the optical properties of the regions are known a priori (Kohlemainen et al. 1999), and later extended to solve simultaneously for both the region boundaries and the optical properties (Kohlemainen et al. 2000a). A further refinement was to use an initial static image reconstruction to identify the number of regions and their approximate locations before finally solving for the boundaries and optical properties of the regions (Kohlemainen et al. 2000b). In 3D, it may be more appropriate to assume that regions of interest are ellipsoidal and to solve for the parameters of the ellipse rather than Fourier components (Kilmer et al. 2003).

A promising technique for the inverse shape estimation problem is provided by representing internal regions using a level set function. This is a function which is 1 inside the regions and 0 elsewhere. This method assumes that the optical properties of the background and regions of interest are known, but does not place constraints on the number of regions or their shapes. Dorn (2002) implemented this method, thresholding an initial spatial reconstruction to provide a starting estimate for the internal objects which are then described using a level set whose boundaries evolve as the reconstruction progresses.