The insertion of genes coding for bioluminescent compounds allows the labelling of cells of interest such as cancerous cells. Such labelling typically provides a low output of bioluminescence photons but is specific, leading to very little background signal. Bioluminescent sources on the surface of a subject can be straightforwardly localised and quantified, but due to the optical scattering and absorption of tissue, the localisation and quantification of an internal bioluminescent source from surface measurements is non-trivial and it is necessary to perform Bioluminescence Tomography (BLT). BLT uses prior information about the optical properties of the subject to reconstruct the location and concentrations of internal bioluminescent sources from surface fluence measurements.

The tomographic reconstruction process requires that the propagation of bioluminescent light through tissue can be modelled. This can be accomplished through analytical and numerical methods [

1]. If the image space (the spatial distribution of bioluminescence in the subject) is discretised into a number of finite and disjoint regions (analogous to pixels in a digital photograph) then the linearity of light propagation allows the propagation model to be formulated as a system of linear equations. This system can be expressed as

**Jx** =

**y**, where

**J** is a matrix representing the physics of light propagation in tissue [

2],

**x** is a vector representing the bioluminescent source distribution (the image), and

**y** is a vector representing the resulting light distribution (the surface measurements). The

**J** matrix is commonly called a Jacobian or weight matrix, and is constructed from prior knowledge of the optical properties of the subject.

Given knowledge of **y** through measurement, and prior knowledge of **J**, it is necessary to reconstruct **x**. However, BLT problems are typically strongly underdetermined as a result of limitations on the number of measurements, and correlation between measurements. For example, in the experimental data examined in this paper it is necessary to reconstruct 10171 unknowns from 526 measurements.

Restrictions on the number of measurements arise primarily from the imaging setup which itself is limited by the duration for which a subject can be sedated, coupled with the typically low signal strength of bioluminescence which necessitates long measurement integration times. Strong correlations between measurements arise from the diffusive nature of light propagation in tissue and the practical constraint of taking measurements only on the surface of the subject. For these reasons, if more measurements are acquired at a particular wavelength, the set of measurements tends to become more correlated, limiting the benefit of increasing the number of measurements. In order to produce accurate images it is then necessary to develop and use reconstruction methods that are robust in the presence of this degree of non-uniqueness. Given small numbers of measurements in comparison to the number of unknowns, prior knowledge about BLT can help reduce the problem of non-uniqueness and so guide the reconstruction process. For example, BLT images the concentrations of bioluminescent chemicals in tissue. As concentrations are non-negative, constraining the unknowns to be non-negative halves the solution space in each dimension. In addition, BLT is commonly used in examinations of cancer growth, many of which grow in a fashion that preserves compactness of the bioluminescent tumour. The very low bioluminescence background in conjunction with the small size and compactness of the bioluminescent tumour result in a spatially sparse bioluminescence distribution, which can be used to further reduce the solution space.

Compressive Sensing (CS) allows the reconstruction of images that are sparse in a known representation given fewer measurements than unknowns [

3], provided certain conditions on the measurements are satisfied [

4,

5]. It has the property of producing exact reconstruction with high probability, through maximisation of reconstruction sparsity whilst minimising measurement error. The measurement conditions require measurements to be:

- independent, so that each measurement provides unique information about the image.
- distributed, so that each measurement provides information about the entirety of the image.

If the set of measurements are sufficiently independent and distributed given the sparsity of the bioluminescence distribution being imaged, then there is a unique solution with maximum sparsity that matches the measurements, and that solution is the desired image. Sparsity maximisation can be used to extract the image from the measurements. In the case of BLT, the restriction of measurements to the subject surface and the diffusive nature of tissue limits the independence of measurements, and light absorption by tissue drastically reduces the sensitivity of measurements to bioluminescent sources deep within tissue. However, the large degree of sparsity possible in BLT images suggests that CS may nonetheless provide benefits over other reconstruction algorithms.

Standard reconstruction techniques attempt to solve

Eq. (1):

where

is the L

^{p} norm. In the case discussed here where the solution is expected to be spatially sparse, CS uses the insight that the most spatially sparse solution satisfying

Eq. (1) is the exact solution with high probability (if the problem satisfies the necessary conditions [

4,

5]). This is equivalent to finding a solution to

that minimises the L

^{1} norm ||

**x**||

_{1}, which induces sparsity in the solution. The L

^{0} norm, which essentially counts the number of non-zero components of

**x**, induces sparsity even more strongly, but is a discontinuous function and so more difficult to minimise. In general, L

^{p} norms for

*p* ≤ 1 can be used to maximise sparsity.

Given the presence of noise in the measurements, minimizing the measurement error exactly is undesirable due to overfitting, and it becomes useful to consider L

^{1} minimization and measurement error minimization as objectives to be optimized jointly, with specified weights on each objective, as in

Eq. (2):

where

*λ* is a term that weights the sparsity constraint against the measurement error. It can be shown [

6,

7] that

Eq. (2) is equivalent to solving

Eq. (3):

where the measurement error threshold,

*ε*, is some function of

*λ*.

A number of CS-inspired reconstruction algorithms have been applied to BLT. Lu et al. [

8] used a differentiable approximation of L

^{1} regularisation. He et al. [

9] used an incomplete variables truncated Conjugate Gradient method, to solve sub-problems involving dynamic subsets of the problem unknowns, in conjunction with L

^{1} regularisation. Cong et al. [

10] used a phase approximation simplification of the Radiative transport equation, permissable region approach and a constrained L

^{1} sparsity minimisation algorithm, with the constraint being a measure of the measurement discrepancy. Yu et al. [

11] used L

^{1} regularisation. Gao et al. [

12,

13] used the full Radiative Transport Equation, with both angularly-averaged and angularly-resolved data, in conjunction with both L

^{1} and Total Variation-based regularisation, to solve for sub-problems on meshes of different resolution and with an adaptive region of interest. The solution they proposed to choose the correct balance of L

^{1} and Total Variation weighting required additional prior knowledge of the bioluminescence distribution. He et al. [

14] used an adaptive region of interest approach with L

^{1} regularisation. Liu et al. [

15] used a dynamic quadratic approximation to the L

^{p} norm, and explored how sensitive the reconstruction process is to the particular regularisation weight used. Zhang et al. [

16] used a primal-dual interior-point approach, which does not need a regularisation weight, and employed an algorithm termination scheme which can be noise dependent. However, multiple termination criteria were used with the same termination threshold, and so its relationship to noise is non-trivial.

One issue not fully addressed in previous works [

8,

9,

11–

15] is that of choosing a regularisation/sparsity weighting, which is the subject of the presented work. Sparsity weighting is equivalent to the imposition of a measurement noise constraint, via the equivalence of

Eq. (2) and

Eq. (3). However, the mapping of sparsity weighting to a desired measurement noise measure (such as sum squared error) or vice versa is non-trivial. Sparsity weighting is also important when using regions of interest, as the sparsity of a given bioluminescence distribution is linked to the size and shape of the region of interest. As the region of interest shrinks around the non-zero parts of the bioluminescence distribution, the distribution necessarily becomes less sparse within the region of interest, and accordingly the sparsity weighting should be changed to account for this.

This work uses the iterative solution of a number of instances of

Eq. (2), with decreasing sparsity weighting, to address this problem (section 2). Firstly, a measurement discrepancy termination threshold is given (based on the instrumental noise statistics). An initial sparsity weighting is generated that heavily favours the sparsity term. A solution to this problem is found, and then used as an initial solution to a new problem with slightly smaller sparsity weighting. This process repeats until the measurement discrepancy termination threshold is satisfied, or the measurement discrepancy term is heavily favoured over the sparsity term. This leads to one of two possible end states:

- The measurement discrepancy threshold is satisfied, and so the result is a local solution of Eq. (3).
- The measurement discrepancy term heavily outweighs the sparsity term, and so the result is a local solution of Eq. (1).

This process of solving a number of sub-problems using a previous solution as a starting point is similar to the approach used by Figueiredo et al. [

6] to improve performance by “warm starting”.

This algorithm was tested against two other algorithms, using simulated and experimental measurements, and its performance quantitatively analysed (section 3).