## Abstract

We describe some modelling and reconstruction methods for optical imaging in the macroscopic and mesoscopic regimes. Beginning with the basic model of radiative transport, we describe the diffusion approximation and its extensions. Some linear and nonlinear problems in diffuse optical imaging are outlined, together with some indications of current trends and future directions.

## 1. Introduction

In this article, we provide a short overview of modelling and image reconstruction in diffuse optical imaging (DOI). The emphasis is on demonstrating the relationships between different techniques and in providing a context for comparison with other more well-established imaging modalities.

Medical imaging methods are generally classified as *indirect imaging* methods in the sense that the sought-for image *x* has to be inferred from data *y* through the inversion of a forward model
1.1In many problems, equation (1.1) is a linear process. In optical imaging, in general, it is a nonlinear process, although for some problems a linear approximation is sufficient. In both the linear and nonlinear cases, the forward model represents the application of a measurement operator to a photon distribution *U*, which itself is the result of a propagator (Green's operator ) that depends on the image *x* and a source term *q*
1.2In emission tomography problems, the source *q* is also part of the inverse problem.

## 2. Modelling in optical tomography

### (a) Physical models of light propagation

For simplicity, we discuss the model for steady-state radiance in this section. The development of time and frequency domain is discussed in §2*b*.

#### (i) The radiative transfer equation

The radiative transfer equation (RTE) is a natural description of light considered as photons. It represents a balanced equation where photons in a constant refractive index medium are propagated along rays , with loss owing to attenuation being balanced by emission, and by photons scattering from other directions. In the absence of scattering, the RTE involves only the transport operator parametrized by the absorption coefficient *μ*_{a}
2.1This is a differential equation in only one variable, with solution
2.2that serves as the basis for the definition of the *ray transform*
2.3In the presence of scattering, and with source terms *q*, equation (2.1) becomes
2.4where *μ*_{s} is the scattering coefficient, *μ*_{tr}=*μ*_{s}+*μ*_{a} is the attenuation coefficient and is the scattering operator, which is local and thus non-propagating.

A series solution for equation (2.4) can be formally written as
2.5known as the *method of successive approximation* [1]. Equation (2.5) is interpreted as successive steps of ray-tracing and local reweighting of rays . A stochastic solver for this series is the Monte Carlo method, which propagates a ray to a random distance drawn from an approximation of the survival probability distribution of unscattered photons, followed by drawing a new ray direction from the probability distribution for the kernel of ; this builds an estimate of *U* photon by photon.

Deterministic solvers seek the direct solution for *U*. The calculation of the ray transform in equation (2.3) can be made computationally efficient, especially with recent advances in graphics processing units, and the inversion of is a purely local operation. However, whichever type of solver is used, the convergence of the series depends strongly on the magnitude of *μ*_{s}, such that the inclusion of more than a few terms becomes expensive. Alternatively, equation (2.5) can be used as the basis of hybrid methods as detailed below.

#### (ii) The *P*_{N} approximation and the diffusion model

If we represent the radiance as a set of spatially dependent variables *ϕ*_{ℓ,m}(*r*) with the angular dependency given in terms of (real) spherical harmonics *Y*_{ℓ,m} of degree ℓ and order *m*
2.6we simplify the structure of , giving rise to (*N*+1)^{2} coupled first-order PDEs in the general case, which we write as
2.7

In the diffusion approximation (DA) framework, the radiance is approximated by zero-order () and first-order () spherical harmonics only, with the result
2.8where *Φ*(*r*) and *J*(*r*) are the photon density and current defined as
2.9and
2.10Inserting the approximation (2.8) into equation (2.4) and following the derivation in Ishimaru [2] and Arridge [3] results in a second-order PDE in the photon density
2.11where *D* is now the *diffusion coefficient* that represents a length scale of an equivalent random walk step. Equation (2.11) and its associated frequency and time domain versions, including the telegraph equation, are the most commonly used in DOI. Another form of coupled second-order PDEs is obtained through the *even-parity RTE*, which generalizes the DA derivation to higher orders of the *P*_{N} approximation. It takes the form
2.12where and and are generalized absorption and diffusion operators, respectively [4].

An example of the difference between diffusion and higher order models is shown in figure 1. The domain is a cube of dimension mm, with a uniform attenuation of and a scattering value of (case 1) and (case 2). An isotropic phase function was used in all cases.

Figure 1(i) shows the photon density for case 1 with a source that is isotropic in the angular variable. Here, the *P*_{1} (DA) and the higher order *P*_{7} solutions produce results that are indistinguishable. This is to be expected as within this regime the RTE is well approximated by the DA. Figure 1(ii) shows the photon density for case 2 with a directed source collimated perpendicular to the boundary. In this regime, one would expect the DA to break down. This figure shows that the *P*_{1} solution is unable to take into account the directional nature of the source and the photon density distribution is diffuse. However, the *P*_{7} solution produces a less-diffuse photon density distribution and the directionality of the source can be seen. In figure 1(iii), we still consider case 2, but now the source is directed at 45^{°} to the boundary and again the *P*_{1} solution is unable to take this into account. The *P*_{7} solution produces a similar photon density distribution as shown in figure 1(ii) but now rotated through 45^{°}.

#### (iii) Hybrid methods

Hybrid methods have been of interest since the early developments in optical tomography. Examples are the Monte Carlo diffusion model [5], the radiosity-diffusion model [6], the transport-diffusion model [7] and the radiative transfer-diffusion model [8].

Another approach was developed in Soloviev & Arridge [9], and shown in figure 2. Here, the series in equation (2.5) is taken to the second term, and the remainder of the series is replaced by the DA with the model equation (2.8). This method allows for the fast modelling of light in regions that are weakly scattering but with localized highly scattering inclusions.

### (b) Time versus frequency

Imaging systems based on the measurement of photon numbers (referred to as ‘direct current’ (DC), or ‘continuous wave’ in the DOI literature) are the simplest to develop. However, they are not always appropriate for a number of reasons that include: the difficulty of calibration, the worse-posedness of the forward model and the non-uniqueness of the inverse model. Therefore, many methods are based on either the frequency modulation of the source, leading to measurements of a complex-valued quantity, or the use of pulsed sources, leading to the measurement of a time-varying signal. The latter may be analysed in the frequency domain through Fourier transform, which combines the robustness of time-domain measurements with the simpler modelling of the frequency domain.

In the probabilistic interpretation, we interpret the time-varying data as a *conditional* PDF in the sense that a photon that has arrived at a given point *r* does so in time interval [*t*,*t*+*δt*] with probability *P*_{r′,t}(*t*|*r*)*δt*. By using the properties of the moment generating function (MGF), we can represent the temporal aspect of the PDF from just a small number of moments, which are computationally easy to calculate [10].

Figure 3 shows an example of the time-varying signal obtained for different source–detector separations across a slab with a heterogeneous distribution of optical parameters. The time-of-flight histogram of transmitted photons at each detector was calculated in two ways: (i) by solving the time-dependent system using a fully implicit finite differencing step in time (see [11] for details) and (ii) by solving for the first three temporal moments and deriving the time-varying solution via the Taylor expansion of the MGF. The comparison is virtually perfect despite the grossly heterogeneous nature of the example which precludes the exact specification of a Green's function. The moment-based method is several hundred times faster.

Using time-domain data allows one to selectively use time gates as data for reconstruction. In figure 4, we compare the reconstructions using early and late time gates. The effect is most noticeable on the scattering reconstruction.

### (c) Series approach to the forward model

Green's operator in equation (1.2) is the inverse of an operator such as those given by the transport, *P*_{N} or diffusion models in equations (2.4), (2.7), or (2.11). If the form of is available for some state *x*_{0}, then we can consider a perturbation in parameters to define a *potential operator*
2.13

Under this notation, a series solution for a perturbed state *x*_{0}+*δx* is given by
2.14where is the identity operator.

Equation (2.14) is a general Born series for forward problems based on PDEs. Other series can be considered as the expansion in equation (2.14) applied to a transform of the data. One widely used transform is the logarithmic transform as given in equation (2.3), which leads to the Rytov series. More generally, the appropriate transform can be based on considerations of the appropriate likelihood term of the forward model in a Bayesian framework.

### (d) Numerical methods

The form of Green's operator is not usually available in arbitrary geometries or for general heterogeneous distributions of optical parameters. Numerical methods are sought in this case.

We take the finite-element method (FEM) as an example. Here, the field *U* and parameters *x* are represented in bases
2.15The basis {*u*_{j}} is chosen to reflect the stability and accuracy of the forward model, whereas the basis {*b*_{i}} is chosen to reflect the expected accuracy of the inverse problem. The PDE is discretized to a matrix K(*x*) parametrized by *x* and with matrix elements (in the Galerkin approximation) and the forward map equation (1.2) becomes the matrix equation
2.16By considering the decomposition
2.17where K^{(i)} is the discretization of the potential operator representing the perturbation *b*_{i}(*r*), we can use the properties of matrix inverses to see that
2.18where U* is the *adjoint field* resulting from the back-propagation of the measurement operator. The term on the right in equation (2.18) forms the matrix elements of the discrete representation of the first term in the Born series; the representation for other series is obtained by transformation as before.

An illustration of the difference between an analytical Green's function model and an FEM model is shown in figure 5. It is apparent that the former is a poor approximation of the case where a boundary is present, so that the use of analytical models of the forward problem induces significant *modelling error*. We return to this issue in §3*f*.

### (e) Reconstruction methods

Image reconstruction methods in optical tomography follow similar strategies to classical tomography problems such as computed tomography and single photon emission computed tomography (SPECT), or to related soft-field imaging methods such as electrical impedance tomography and microwave tomography. A taxonomy of algorithms is described by Arridge & Schotland [12], and some comparisons on a test case are presented by Schweiger & Arridge [13]. In the remainder, we focus only on the *regularized output least-squares method*. Here, the reconstruction is treated as the optimization of a model-fitting problem
2.19where is the sum of a negative log-likelihood term for a multivariate Gaussian noise model with covariance Γ_{e} and a negative log prior term *Ψ*(*x*). For nonlinear forward operators, or for non-Gaussian prior densities, equation (2.19) needs to be solved iteratively through the normal equations
2.20where H_{n} is the Hessian of *Ψ*(*x*_{n}) and *A*′,*A*^{′*} are the forward and adjoint Fréchet derivatives, given by the representation in equation (2.18).

For large-scale problems, the inversion step in equation (2.20) is a limiting step. Instead, the Krylov method constructs only a few terms in the series
2.21where {*v*_{k}} are the Krylov basis vectors, and is a smoothing operator representing the correlation of the prior. Equation (2.21) involves forward and backward projections and filtering in the data and solution spaces. It can be implemented in a *matrix-free* framework without the need to construct the memory-intensive derivative matrices *A*′.

## 3. Inverse problems in optical tomography

### (a) Parameter identification problems

The initial formulation of optical tomography was in terms of the recovery of intrinsic optical parameters from boundary measurements of transmitted and/or reflected light. This has had applications principally in brain imaging and in breast cancer detection. It needs to be considered as a nonlinear inverse problem, although for small changes in state, a linearized approximation has also been used to some success. An illustrative example is shown in figure 6, which was solved using the method of §2*e*.

### (b) Linear source identification problems

Several types of imaging problems are associated with source identification. The spontaneous emission of optical radiation from within tissues through the action of luciferase on its substrate, luciferin, leads to bioluminescence tomography (BLT). It is the optical analogue of SPECT in nuclear medicine.

In the study of Gu *et al.* [14], the bioluminescent source distribution was recovered from a monochromatic dataset and theoretical aspects of uniqueness of solutions were discussed in Wang *et al.* [15] and Cong *et al.* [16]. The idea of using multiple wavelengths in bioluminescence imaging was introduced by Chaudhari *et al.* [17]. Combination of bioluminescence and positron emission tomography is discussed in Alexandraxis *et al.* [18].

In fluoresence optical tomography, the emitted light is the result of the promotion of endogenous or exogenous *fluorophores* to an excited state by an incident ‘probe’ field, followed by exponential decay through the emission of fluorescent radiation. The detected radiation is at a longer wavelength (i.e. lower energy) than the radiation used as the source. Because the excitation is under the control of the experimenter, and can be constant, frequency or time domain, a richer set of data is available than in BLT.

The inverse problem is concerned with the recovery of the density and/or lifetime of the fluorophores. The latter requires time or frequency domain measurements. In the frequency domain, we represent the parameter of interest as
3.1where *η*_{0} is concentration and *τ* is lifetime. In the case where the background optical properties are known for both the excitation and fluorescent wavelengths, the inverse problem is linear, as it is for BLT.

We show in figure 7 a target and its BLT reconstruction. The internal fields are shown for six different wavelengths, where the increasing absorption leads to more rapid damping of the photon density as a function of distance from the target. However, these fields are highly correlated which emphasizes the strong illposedness of this problem. By contrast, in figure 8, the same target is considered as a fluorophore and its reconstruction can be seen to have a much higher resolution. In this case, the internal fields are the product of the stimulation field, controlled by the experimenter, and the fluorophore concentration; as can be seen, the correlation of these fields is much less, which leads to a less severely illposed inverse problem. In both cases, the reconstructions were made with the matrix-free Krylov method with Perona–Malik edge-preserving regularization.

### (c) Multi-spectral methods

Although each of the above problems could be considered at a range of spectral samples and the spectral variation of the recovered solutions determined, the idea in multi-spectral DOI is to reformulate the problem into the recovery of a set of images of known chromophores with well-characterized spectral dependence:
3.2where ** ϵ** is a known matrix. One phenomenological model of the spectral dependence of scattering, based on an approximation to Mie scattering that has proved useful, is
3.3We note that the above model for the wavelength dependence of

*μ*

_{s}′ corresponds to Rayleigh scattering when

*b*=4. In general, subdominant power-law corrections may be necessary to accurately represent the scattering behaviour of tissue.

The use of the spectral constraints (equations (3.2) and (3.3)) allows for a reparametrization from a set of uncoupled inverse problems (one for each spectral sample) into a problem in recovery of the parameters {*a*,*b*,*c*_{k};*k*=1,…,*K*}. As well as reducing the inverse problem dimension, this approach can ameliorate some of the aspects of non-uniqueness between absorption and scattering in the case of DC measurements, because the necessary conditions for functions to be in the problem null-space are invalidated under the spectral constraint conditions [19,20].

Figure 9 shows multi-spectral reconstructions from simulated frequency domain data at two wavelengths in a circular domain. Scattering parameter *b* was set to a homogeneous target value of *b*=1, while an inhomogeneous target distribution was used for *a*. The target distributions for the two chromophores and scattering parameters are shown in figure 9*a*. The resulting absorption and scattering coefficient distributions at the two wavelengths are shown in figure 9*b*. The reconstruction results are shown in figure 9*c*.

### (d) Choice of priors

A general form of prior is defined by the functional
3.4Here, the function *ψ* controls the penalty associated with edges in the reconstructed image *x*, weighted with respect to an orientation metric W reflecting an associated image, assumed to be well defined. In a multi-modality framework, W is an edge indicator which is chosen to enforce smoothing tangentially to edges of the associated image, and to decorrelate smoothing across them. A simple example is shown in figure 10.

### (e) Data compression and structured light

In optical imaging applications for small targets, such as mice, the amount of transmitted or emitted photons is large and the use of cameras as detectors has become prevalent. This brings with it the problem of a large dataset, which precludes the explicit storage and inverse of the forward operators or their derivatives. The use of matrix-free methods removes the memory-intensive aspects, but still can lead to slow inversion. This has led to the idea of the projection of the recorded camera data into a compressing representation, such as a Fourier or wavelet basis. Similarly, it is natural to consider distributed illumination of the input surface, and a compressive basis for the source; this has come to be known as the *structured light* method.

In the particular case of an infinite slab and the use of Fourier patterns, the structured light technique corresponds to the analytical block-decomposition method of Markel & Schotland [21], which has been successfully applied to the reconstruction of complex structures in measured data [22,23]. An illustration is shown in figure 11, where an FEM model of a slab is used as the object.

More generally, light patterns can be projected onto, and captured from, an arbitrary geometry using numerical methods familiar in computer graphics. First results in a cylinder using wavelet sources and patterns have recently been reported [24].

### (f) Approximation error modelling

In the approximation error method [25], we abandon the need to produce an ‘exact’ model. Instead, we attempt to determine the statistical properties of the modelling errors and compensate for them by incorporating them into the image reconstruction using a Bayesian approach. The effect of the modelling error can in many cases be of significantly higher order than that of the measurement noise. By characterizing its statistics, we can compensate for bias and nuisance correlation in the images to a remarkable degree. In the study of Arridge *et al.* [26], this technique was shown to result in reconstructed images using a relatively inaccurate forward model that were of almost equal quality to those using a more accurate forward model; the increase in computational efficiency was an order of magnitude.

We show a different application in figure 12. Data were generated in an accurate three-dimensional FEM model but reconstructed using a linear model constructed from analytical Green's functions [27]. Subtraction of the data from a homogeneous model provides a simple non-Bayesian model error correction. We see that the use of an analytical model for inversion fails even when using a reference measurement to calibrate data, if the latter is taken under incorrect estimates of the reference optical properties. The use of the approximation error method compensates for the incorrect reference estimate.

## 4. Conclusions

In this article, we gave an overview of strategies and recent advances for DOI. Issues in both modelling and reconstruction were considered. We limited the discussion only to optical imaging, mainly in the diffuse regime, with some indication of the extension into the weakly scattering regime, sometimes known as ‘mesoscopic imaging’. Space precludes several interesting topics such as techniques exploiting coherence or polarization. Currently, one of the most active areas is the use of ‘coupled physics’, notably the interaction of optical and acoustic waves. We believe that several of the techniques for optical image reconstruction can be of use in these applications, in particular, for the quantitative recovery of optical parameters from recovered acoustic source images.

## Acknowledgements

The author thanks T. Correria, S. Prerapa, M. Schweiger, V. Soloviev and T. Tarvainen for providing several of the illustrations. Funding was provided by EPSRC grant EP/E034950/1. Modelling and reconstruction results were generated using the TOAST package (http://www.toastplusplus.org).

## Footnotes

One contribution of 20 to a Theo Murphy Meeting Issue ‘Illuminating the future of biomedical optics’.

- This journal is © 2011 The Royal Society