Gradient-based quantitative image reconstruction in ultrasound-modulated optical tomography: first harmonic measurement type in a linearised diffusion formulation

Ultrasound-modulated optical tomography is an emerging biomedical imaging modality which uses the spatially localised acoustically-driven modulation of coherent light as a probe of the structure and optical properties of biological tissues. In this work we begin by providing an overview of forward modelling methods, before deriving a linearised diffusion-style model which calculates the first-harmonic modulated flux measured on the boundary of a given domain. We derive and examine the correlation measurement density functions of the model which describe the sensitivity of the modality to perturbations in the optical parameters of interest. Finally, we employ said functions in the development of an adjoint-assisted gradient based image reconstruction method, which ameliorates the computational burden and memory requirements of a traditional Newton-based optimisation approach. We validate our work by performing reconstructions of optical absorption and scattering in two- and three-dimensions using simulated measurements with 1% proportional Gaussian noise, and demonstrate the successful recovery of the parameters to within +/-5% of their true values when the resolution of the ultrasound raster probing the domain is sufficient to delineate perturbing inclusions.


Introduction
The wavelength-dependent optical absorption and scattering coefficients of biological tissues provide clinically valuable information regarding tissue function and composition. Purely optical techniques such as diffuse optical tomography (DOT) are capable of measuring these coefficients but suffer from limited spatial resolution due to the high degree of optical scattering encountered in typical biological media [1][2][3]. Ultrasoundmodulated optical tomography (UOT) is a hybrid technique which aims to recover the coefficients with significantly improved resolution by combining the optical contrast of near infra-red light with the spatial resolution of focused or timegated ultrasound fields.
Much effort has been expended in advancing the experimental technique in UOT. The problem of detecting the small ultrasound-modulated optical flux against the large unmodulated background has received significant attention. This problem is particularly challenging since the requisite use of a coherent source generates a spatially incoherent speckle pattern on the boundary of the domain. The flux of individual coherence areas must therefore be collected in parallel [4][5][6], or manipulated such that their contributions can be measured in summation [7][8][9][10][11].
Less attention has been paid to the fundamental problem that hybrid techniques such as UOT are only capable of producing quantitative images under some form of model-based reconstruction procedure. This was succinctly demonstrated by the images produced by Lev and Sfez [12,13] where the raw data from a UOT measurement can be seen to resemble the optical sensitivity functions [14] in a given domain. UOT, unlike alternative hybrid techniques such as photo-acoustics (PAT), has similar sensitivity to both absorption and scattering perturbations [15]. Raw measurements attempting to quantify absorption perturbations will be significantly affected if an erroneous assumption of constant optical scattering is made.
Previous investigations have successfully recovered the optical absorption coefficient from simulated [16,17] and experimental data [18]. Simultaneous recovery of both the optical absorption and scattering coefficients raises the more subtle problem of non-uniqueness. The key idea is that the recovery of the two coefficients requires at least two sets of internal data, as has been demonstrated in a related incoherent formulation of UOT [19,20]. We recently demonstrated that uniqueness can be restored by the use of multiple optical source and detector locations [21].

Overview and contribution
This work is organised as follows. In section 2 we provide an overview of forward modelling techniques for UOT, and derive an efficient model of the power-spectral density of the UOT signal from a non-linear time-domain form presented elsewhere in the literature. In section 3 we pose the inverse problem of recovering the internal optical parameter distributions from measured data, before implementing our adjoint-assisted gradient based optimisation technique in section 4. We employ this optimisation in section 5 by performing, for the first time, simultaneous reconstruction of the optical properties in a simulated UOT experiment in two-and three-dimensions. We close this work with a discussion of our findings in section 6.

The physical basis of UOT
As an acoustic wave propagates through a biological medium it induces small changes in the refractive index of the medium [22], and causes the displacement of optical scatterers from their rest position [23,24]. Under coherent illumination an otherwise static (in the absence of Brownian motion) speckle pattern is generated on the surface of the medium which changes in time as the optical path lengths of the scattered waves travelling through the medium are phase modulated by the acoustic field. These changes may be measured as either a temporal decorrelation of the intensity autocorrelation function, or, equivalently, as the modification of the power spectral density of the measured light.

Forward modelling techniques
Various models of this process have been presented in the literature. The starting point in each case is a time-domain description of the phase perturbations applied to the scattered optical waves propagating through the medium.
Path-integral methods develop expressions for the total average phase perturbation over an optical path of a given length by averaging the phase perturbation expression over all free-paths and scattering events. Under an assumption of weak scattering the perturbations applied to each optical path length are considered independent, and each provides an individual contribution to the total field autocorrelation function. Integration of the contributions to the correlation function over a probability distribution of path lengths, typically found analytically from the diffusion equation, yields an estimate of the measured field autocorrelation function. This approach is similar to some of the original investigations into diffusing wave spectroscopy by Maret and Wolf [25], and Pine et al. [26].
The principal limitation of path integral methods is that they can only incorporate planar acoustic fields (owing to the averaging over all potential paths). Despite this limitation, the technique has offered significant insight into UOT for time-harmonic acoustic fields in isotropically [22] and anisotropically scattering [27] media, and for acoustic pulses [28] in which there are significant correlations in the phase perturbations between successive scattering events.
Correlation transport is an extension of radiative transport theory (a high-frequency approximation for optics) to consider media in which there is a temporal variation of the underlying medium. This idea was originally investigated by Ackerson et al. [29] who proposed a correlation transport equation (CTE) for use in the field of diffuse correlation spectroscopy (DCS). Dougherty et al. [30] later provided a more rigorous derivation based upon analytic theory with moving scatterers, building on the work of Ishimaru and Hong [31] and others [32,33]. A similar approach can be taken in UOT [34,35] where the scattering objects now move deterministically, and additional terms are introduced to account for the modulation of the refractive index.
The only necessary assumption in a CTE based model is that of weak scattering (l tr λ, where l tr = (µ a + µ s ) −1 is the optical mean free path, µ a is the optical absoprtion coefficient, µ s = µ s (1 − g) is the reduced scattering coefficient, g is the scattering anisotropy, and λ is the optical wavelength). Owing to the complexity and high-dimensionality of the phase space of the integrodifferential form of the CTE in UOT, solutions have to date only been presented using statistical approximations sought via the computationally expensive Monte-Carlo (MC) method. We have previously presented GPU implementations of a CTE for UOT [36,37] which demonstrate significant improvements in speed over traditional implementations, but even so, the computation requirements are such that at present these models are only suitable for use as a 'gold-standard' by which more approximate techniques may be validated.
Correlation diffusion is an approximation of correlation transport by a first order expansion of a given CTE in spherical harmonics, valid under a number of assumptions, principally that µ s µ a . The diffusion approximation is readily made for the RTE in the context of DOT, and also the CTE which arises in DCS, see Boas [38] for a full derivation. The process is more complicated in UOT, principally due to correlations between phase increments at successive scattering sites. By combining aspects of the statistical averaging of the path-integral methods with the transport component of the CTE, Sakadžić and Wang [39] developed a correlation diffusion equation (CDE) for the field autocorrelation function in UOT.
Approximations in the derivation of the CDE presented in [39] require conditions of weak modulation (limiting the acoustic pressure amplitude to circa 10 5 Pa for typical medical ultrasound frequencies), and uncorrelated phase increments (k a l tr 1, where k a is the acoustic wavenumber), in addition to the assumption of weak scattering inherent to the underlying CTE. The flexibility and computational simplicity of this model are highly attractive for application in an image reconstruction application.

Linearised power-spectral correlation diffusion model
The CDE presented in ref. [39] describes the wide-sense stationary optical field autocorrelation function φ(r, τ) in a given domain, and it can be written: with a lossy diffusion operator where κ = (3µ s ) −1 , is the diffusion coefficient, q 0 (r) is a coherent isotropic source source, A is related to the index of refraction mismatch between the turbid and external media [40], n is a unit vector normal to the boundary of the domain, ω a is the acoustic angular frequency, τ is lag, η(r) ∝ P 0 (r) 2 is the acousto-optic modulation efficiency, which has weak dependence upon the optical parameters, and P 0 (r) is the local acoustic pressure amplitude. In previous work [17] we employed this approximation as part of a linear reconstruction technique for UOT, neglecting the weak dependence upon optical coefficients in η(r), successfully recovering absorption coefficients from boundary measurements. This model is readily applicable to techniques in which the intensity or field autocorrelation function are directly recorded in the temporal lag domain. However, many interferometric detection methods in UOT directly record the powerspectral density of outgoing modulated flux at the acoustic frequency (the first harmonic flux). To calculate this value multiple runs of the forward model must be made at difference values of lag, τ. Given that the magnitude of η(r) is already limited by the assumption of weak modulation, we are motivated to further linearise the model of equation 1 to arrive at a direct frequency (power-spectral) domain representation.
Suppose we take two measurements of the same medium under two acoustic pressures such that the field correlation function φ(r, τ) → φ b (r, τ) + φ δ (r, τ), under η → η b + η δ , then Subtracting from equation 4 the expression for the 'baseline' measurement made under η b results in a non-linear expression for the perturbed field under changing insonification, We denote by R(η) the forward operator which returns φ(r, τ) in a given domain by implementing equations 1 and 2 under parameterisation by η. We employ this notation to rewrite equation 5, where we identify the first Fréchet derivative of the non-linear forward operator, R , and the second order terms in the small parameter η δ which we neglect to complete our linearisation: This form has significant value since judicious choices of ultrasound pressures and a known pressure squared dependence may permit η(r) to be identified experimentally. For now, we choose our background measurement to be made in the absence of an acoustic field such that η b = 0, and η δ = η. Equation 6 becomes R(η) − R(0) + O(||η 2 δ ||) = η δ R (0), and after dropping the second order terms we isolate the coupled expressions for the unperturbed field, and the perturbation, By inspection, it is evident that φ(r, τ) = φ b (r, τ) + φ δ (r, τ) contains spectral content only at DC and the frequency of the ultrasonic excitation ω a . We may now find an expression for the power-spectral density of the fluence in the domain by the Wiener-Khinchin theorem. Taking the Fourier transform of the field autocorrelation function we find, where we take φ i (r) to mean φ(r) or φ 1 (r), φ 1 (r) =φ(r, ω a ) is the first harmonic correlation fluence, and the total fluence rate, equivalent to the CW fluence in a DOT experiment We also define that the fundamental (DC) correlation fluence φ 0 (r) = φ(r) − φ 1 (r) =φ(r, 0). The information contained in φ(r) has limited spatial localisation due to location of the source term q 0 (r). The components at the fundamental and first harmonic, φ 0 (r) and φ 1 (r), both contain spatially localised information related to the total fluence rate in the medium by the presence of the modulation efficiency term. In this model we may readily identify the oft-cited 'virtual acousto-optic source' as the product η(r)φ(r). We define our measurement as the power spectral flux across the boundary at the ultrasound frequency. For a point detector this is given by [2], If we take the biological medium under investigation to be index matched to its surroundings, then the constant A in the Robin boundary condition takes a value of unity. Combining the expression of the flux with the boundary conditions gives a simple form for the measurement: We introduce a measurement operator, M, which integrates the flux over a finite region of the boundary of the domain, a projection operator, P, which incorporates this process, and define an operator F which implements equation 11 under parametrisation by the optical parameters of interest, Before proceeding to consider the inverse problem, we make two comments on our model. First, we note that if η(r) = δ(r − r 0 ), the 'internal data' returned by a UOT measurement is in fact the fluence in the domain at r 0 propagated by the Green's function which solves the underlying diffusion approximation. This suggests that some form of 'inverse diffusion' approach to the image reconstruction problem may be fruitful. We will return to this point in the discussion. Second, a model similar to that of equations 8 and 9 was previously developed by Allmaras and Bangerth [16] by the application of the Born approximation to the modulated field in a path integral formulation, the nature of those paths then being formalised in the diffusion framework (similar to earlier works by [27]). Our derivation shows that this model can be derived by reasoned approximations to a correlation transport equation, and that it is then amenable to a power-spectral representation.

Measurement protocol and notation
To this point we have considered the case of a single optical source and detector, and ultrasound field distribution. To perform imaging we will probe the domain of interest with multiple insonification profiles and optical source-detector pairs. We define a measurement index For a given experiment the data vector y ∈ R N m , N m = N i × N j × N k contains data from all combinations of sources, detectors and acoustic field profiles, and this is found by the previously introduced continuous to discrete projection and measurement operators of equation 17.
To refer to a single measurement the subscript ρ implies a particular choice of the 3-tuple {i, j, k} which corresponds to the appropriate subset of the forward operators, and potentially their derivatives. That is to say that, For convenience we also denote the compound optical parameters,

The inverse problem
To perform a reconstruction we begin by making the assumption that our measurements y m = y + n are corrupted by noise n ∼ N(0, Γ e ) drawn from a normal distribution with mean zero and covariance Γ e . Given a known parameter distribution, the probability of the data is thus given by We assume that the parameters we seek to reconstruct vary smoothly in space such that Lx ∆ 2 is small, where (x ∆ = x − x g ), x g is a suitable a priori reference parameter set, L =x −1 g ∇, andx g = (μ a,g ,μ s,g ) T is the mean of the compound reference parameter distribution. This corresponds to the probability distribution of the parameters [41] π with the improper covariance Γ x =x 2 g ∇ −2 . This corresponds to a classical choice of first order Tikhonov regularisation, which promotes smooth solutions, where the termx g serves to compensate for the different scaling between the two parameters, sphering the solution space.
From a Bayesian perspective, the image reconstruction problem corresponds to finding some appropriate metric of the a posteriori probability distribution π(y|x)π(x).

The error functional, and its gradient
In this work we define the parameters of our reconstructed image x * as the maximum a posteriori (MAP) estimate of π(y|x)π(x). This corresponds to a minimisation of the weighted least squares error functional [1] x * = argmin where we have introduced an ad-hoc hyper-parameter λ which serves to control the relative contribution of the data term and the regularisation term to the error functional. Under the assumption that the error functional is convex, minimisation of E(x) corresponds to the solution of the nonlinear equation E (x) = 0. To find this expression we begin by expanding equation 24, denoting the residual b = (y m − P[x]) and taking the derivative with respect to x, where P [x] is the Fréchet derivative of the projection operator, The functions C[x] define the sensitivity of the measurements to perturbations in their respective optical parameters. Consistently with our previous work [17], we refer to these functions as the fundamental correlation measurement density functions (CMDFs) in analogy with the associated photon measurement density functions arising in DOT [14,42].

Implementation
To proceed we must now redefine our problem in a finitedimensional setting suitable for solution using numerical methods, and decide upon a technique by which we solve the large non-linear system resulting from the discretisation of the expression E (x) = 0.

Forward model
The domain under consideration is subdivided into a mesh of non-overlapping elements joined at N n vertex nodes. On this mesh we define a set of piecewise linear basis functions such that u i (r j ) = δ i j for i, j = 1, . . . , N n where r j located at the j th vertex node. We may then define finite dimensional approximations to the parameter distributions, solutions, and sources, where we take χ k to represent the k th component of the N n × 1 vectors χ of nodal coefficients which define the discrete approximationχ(r) of any of the continuous functions µ a , µ s , κ, η, φ 0 , φ 1 , q 0 , q 1 , and x = (µ a , µ s ) T is the 2N n × 1 compound vector of optical coefficients. We solve the correlation diffusion equation by the finite element method [42][43][44]. Equation 11 is multiplied by test functions which obey the boundary conditions, and whose zeroth and first derivatives are integrable over the domain. The boundary conditions of equation 13 are incorporated by subsequent integration by parts. Selecting the test functions in the weak formulation to be the same as the basis in which we have defined our parameters allows us to write the resulting linear system, are the N n × N n matrices of system matrix integrals, formed by the basis system matrices and To find the vector of virtual source coefficients q 1 we begin by identifying the discrete form of the right hand side of equation 11: which we insert into equation 34, The vector of coefficients φ employed in the definition of q 1 themselves come from the solution of the system which is the standard formulation of the continuous-wave DOT problem with source term coefficients found by projection of the source profile into the basis. In this work we employed the Toast ++ toolbox to manage mesh data, peform the relevant elemental integrals, and assemble the system matrices and righthand-sides [45].

Error functional
The discrete form of the error functional in equation 25 is found by replacing the parameters by their nodal approximations, and the non-linear projection operator P with the discrete form found via the finite element method, dr, (38) where for a given measurement Before we derive the form of the gradient of the error functional, we will find the matrix form of the regularisation term.
Taking µ to represent either µ a or µ s , with µ ∆ = µ − µ g , such that, where Thus, regularisation applied to the compound vector x = (µ a , µ s ) T can be written, where

Error functional gradient
We may now take the derivative of the discrete error functional of equation 38 with respect to each of the parameter coefficients, and insert our definition of the residual, The derivative of the projection operator P [x] = J[x] is the dense N m × 2N n Jacobian matrix of partial derivatives. This matrix, which constitutes the discrete approximation to equation 27, consists of pairs of discrete correlation measurement density vectors C µ a ρ , and C µ s ρ for each parameter µ a and µ s , for every measurement ρ. We take the derivative of equation 39 for measurement ρ with respect to the i th nodal value of the parameter µ, where we understand the measurement operator, fluences and sources to be those corresponding to measurement index ρ. The derivatives of the system matrix were defined in equations 31 and 32. We proceed by taking the transpose of equation 46, before defining φ + as the solution to the adjoint equation Sφ + = M and expanding the transposed term, Recalling the definition of q 1 and noting the symmetry of the basis system matrices, we may take the derivative, which, assuming the coherent optical source q 0 (r) which generates the field φ is independent of the optical properties, can be written, We now define the adjoint first harmonic source which is propagated to form the adjoint first harmonic fluence, Taking equation 50 in product with φ + , and inserting our definition of φ + 1 , leads to a final form of a given CMDF We note that the form of the sensitivity function corresponds to an approximate adjoint method we previously proposed for use in MC simulations [37]. Inserting equation 54 into equation 45 for all measurements and nodal parameters gives an expression for the functional error gradient, By employing the adjoint solutions to the forward problem we have derived a method by which the Jacobian can be calculated with at most four runs of the forward model for each measurement index, ρ. In a UOT experiment, it is highly unlikely that each ρ will consist of an entirely unique set of optical sources, detectors, and acoustic field distributions. In this case pre-computation of all optical solutions and adjoints φ and φ + will result in significantly accelerated calculations.
A typical approach to minimising E(x) is via a Newton based method such the Quasi-Newton, Gauss-Newton and Levenberg-Marquardt techniques. Whilst such approaches avoid computation of the Hessian matrix by various approximations, the underlying Taylor expansion of the objective function still calls for the storage (and repeated inversion) of the complete Jacobian matrix. In techniques such as DOT there are typically a limited number of sources and detectors such that the viability of storing the Jacobian is principally dependent upon the number of degrees of freedom in the forward model, as determined by the finesse and dimension of the discretisation: three-dimensional problems typically have a far greater number of degrees of freedom (DOF) than their two-dimensional counterparts. This is also true of UOT, however in this application there may be an arbitrary number of measurements corresponding to different acoustic field profiles such that Jacobian matrix extends in both the row-and column-space. We therefore take an alternative approach which is to directly calculate the functional gradient on a row-by-row basis, and employ this in a gradient based optimisation technique which does not require storage of the explicit Jacobian. This technique has been previously demonstrated in DOT [46] and quantitative photoacoustics [47,48].
To proceed we incorporate the error covariance matrix into the residual such that b * = Γ −1 e b, and rewrite equation 56, We define two adjoint equations which correspond to those conceived in the derivation of the Jacobian matrix, with modified adjoint first harmonic source term and the N n × 1 vector z u Equation 61 permits the direct calculation of the gradient of the error functional without explicitly building the intermediary Jacobian matrix.

Optimisation
To solve the non-linear system E (x) = 0 we employed a (Polak-Ribiére) non-linear conjugate gradient method. The gradients were preconditioned with the block mass matrix: For details of this optimisation technique, refer to standard texts, e.g., Nocedal and Wright.

Results
In this section we will examine the form of the CMDFs defining the sensitivity of our measurements, and demonstrate simulated reconstructions in two-and three-dimensions.

Correlation measurement density functions
The CMDFs which define the sensitivity of our measurement to perturbations in the optical properties of the system offer significant insight into this imaging modality. We will now consider the form of the individual CMDFs, under variation of the acoustic field and optical source-detector profiles.
For this purpose we consider a two-dimensional circular domain of diameter 50mm with homogeneous optical properties. The absorption and reduced scattering coefficients µ a = 0.01mm −1 , µ s = 1mm −1 , of the domain are typical of biological tissues. The refractive index of the domain is matched to its surroundings. Around the periphery of the domain are placed three optical sources and three optical detectors two of which are collocated, each having Gaussian profiles of full-width halfmaximum (FWHM) 5mm. Given the reciprocity of the problem their exist a set of six source-detector combinations which return unique data. The domain is discretised into a set of 6, 840 linear triangular elements joined at 3, 511 vertex nodes. A set of 349 focused acoustic fields probes the domain through the two-dimensional plane. Each field has a Gaussian profile with FWHM = 2mm, and the set are arranged over a rectangular grid with spacing 4mm, truncated at a radius of 22mm from the centre of the domain. The peak magnitude of η(r) = 0.25. The discretised domain is depicted in figure 1, where each circle represents an acoustic focal point.
An expression for the discrete form CMDFs was presented in equation 54. This expression consists of four terms, 1. φ, the forward total fluence, equivalent to the CW DOT fluence in the medium due to application of a given optical source. 2. φ 1 , the first-harmonic UOT fluence which results from the virtual acousto-optic source given by the product of the light distribution from the optical source and the acoustooptic efficiency term. 3. φ + , the adjoint total fluence, equivalent to the CW DOT fluence in the medium due to an adjoint source found by application of a given optical measurement operator. 4. φ + 1 , the adjoint first-harmonic UOT fluence which results from the adjoint virtual acousto-optic source given by the product of the adjoint fluence distribution from the adjoint optical source and the acousto-optic efficiency term.
In figure 2 we plot each of these terms for a specific acoustic field distribution and optical source-detector pair. The sensitivity function in µ a consists of the summation of two products, which we also depict. The first product φ × φ + 1 (depicted in the fifth image of figure 2) represents the sensitivity of the measurement due to attenuation of the input illumination which generates the first harmonic modulated fluence. We see a peak sensitivity in the region of the acoustic focus: this is partially due to our choice of a distributed optical source. Evidently if the optical source were point-like, the sensitivity to a perturbation in µ a near this point could potentially be larger than that in the region of the distributed acoustic field, since a reduction in the strength of the optical source will cause a linear change in the total fluence. The second product φ 1 × φ + (depicted in the sixth image of figure 54) represents the sensitivity of measurement due to attenuation of the modulated field prior to detection. As before, the exact form of this term is determined by our choice of detector profile. In summation, these two terms represent the sensitivity of a measurement of the modulated fluence due to a perturbation in the optical absorption profile. These figures   demonstrate significant similarities to those derived for the nonlinearised form of the forward model which we investigated in [37]. However the use of a point source and small aperture detector in the referenced work lead to increased sensitivity near the source and detector regions (we will return to this point in the discussion).
In figure 3 we plot the CMDFs for µ a and µ s which arise for all six source detector pairs and a given ultrasound field distribution. The form of the sensitivity functions of the absorption coefficient follows our previous exposition; in each case a region of peak sensitivity is seen near the acoustic focus, extending outwards towards the given source and detector position for that measurement. The form of the sensitivity functions for µ s retain the obvious dependence on the source and detector location, but have a more complicated structure in the region of the acoustic focus by virtue of their dependence upon the divergence of the four fields from which they are generated. In particular, we see a (spatially) fast reduction in sensitivity in the region where the forward and adjoint optical sensitivities are shadowed by the acoustic focus. We find areas of negative sensitivity in the shadowed region, where an increase in scattering will cause more of the input light to be modulated and detected.

Two-dimensional reconstruction
To demonstrate our algorithm in the reconstruction of an image we introduce perturbations in the absorption and scattering coefficients of the two-dimensional domain utilised in section 5.1. Simulated measurements were performed for all sourcedetector pairs and 1% proportional Gaussian noise added to the data. The error covariance matrix Γ −1 e was set to diag(1/y) in accordance with our noise model. The regularisation hyperparameter was selected by inspection. The optimisation was performed by the non-linear conjugate gradient method, and was terminated when the change in the objective function fell below ∆E (x) ≤ 0.01, which resulted in convergence after 26 iterations. Figure 4 depicts the target optical parameters, the results of our reconstruction, and the percentage error for each coefficient. The reconstructed images can be seen to be in excellent agreement with their associated targets. The error in µ a appears to be correlated with the image insofar as there is a slight under-reporting of the peak-positive and negative perturbations. The error does not exceed ±5% at any point in the image. A similar result is seen in the error in the reconstruction of µ s .

Three-dimensional reconstruction
To demonstrate a three-dimensional reconstruction we consider a cylindrical domain of diameter and height 50mm. The absorption and scattering coefficients consist of an homogeneous background µ a = 0.01mm −1 , µ s = 1mm −1 , with numerous perturbations in both parameters with a maximum magnitude of 50% of the background value. A set of 1456 focused acoustic fields with three-dimensional Gaussian profiles of FWHM = 3.5mm probe the domain. The focal points are arranged over a rectangular grid, truncated at a radius of 22mm.
The peak magnitude of η(r) = 0.25. The geometry of the problem, the location of the focal points of the ultrasound fields, and isosurfaces of the optical parameter perturbations are depicted in figure 5. Four sources and detectors were arranged around the periphery of the domain. The first three sources and detectors were located as per the two-dimensional problem, at z = 0. The final source and detector were placed on the top and bottom surfaces of the domain, respectively. Each source had a profile corresponding to a cosine window of diameter 20mm. Simulated measurements were performed for all source detector pairs and 1% proportional Gaussian noise was added to the data. The error covariance matrix and regularisation hyper-parameter were determined as per the two-dimensional case. The optimisation was performed using the preconditioned non-linear conjugate gradient method, and was terminated when the change in the objective function fell below ∆E (x) ≤ 0.01, which resulted in convergence after 61 iterations. Figure 6 compares the target optical parameters, the results of the reconstruction, and the percentage error for each coefficient for a slice through the domain.
The reconstructed images successfully capture the location and magnitude of the larger perturbations in the optical coefficients of the domain with excellent accuracy. The smaller perturbations are somewhat under-reported, especially in the case of the absorption coefficient.
This can be directly attributed to the use of a coarser grid of ultrasound focal points, and the increased FWHM of the fields, which in this case is comparable to the size of the smaller inclusions. As such, this system is significantly more underdetermined than the two-dimensional case, as the number of degrees of freedom has grown significantly compared to the number of measurements (ultrasound field locations).  The result of this is that the weighting of the prior, which in this case enforces smoothness, is relatively larger. This leads to the small absorption perturbation 'merging' into the larger feature, which can be seen both in the isosurface and slice plots. An edge preserving prior such as total variation may better retain delineation between regions of differing optical properties.

Discussion & Conclusions
Before closing, we consider some points raised by this work, and make suggestions for future investigations.

The coherent optical source, and uniqueness
We noted earlier that the use of a diffuse source term of similar spatial extent to the detector profile resulted in UOT sensitivity functions with a maximum sensitivity at the point of the acoustic field, and that this differed from previous work [17] in which a point source dominated the sensitivity of the system for this measurement type. In another recent publication which employed point sources [21] we noted that interchanging of the source and detector locations lead to an improvement in the spectrum (the number of significant singular values) of the approximated Hessian which arises in a quasi-Newton optimisation: this has significant implications regarding the uniqueness of the reconstruction. This previously unexplained phenomenon seems at odds with the symmetry which results from the physical reciprocity of the system.
In this work, with identical source and measurement apertures, no advantages were gained from taking measurements with swapped source-detector locations. Thus, the extra information found by transposing the source and detectors in [21] was in fact achieved by virtue of their differing profiles. Whilst this might initially suggest an advantage in using a point source and diffuse detector, we must consider that the extreme sensitivity near the source position may be deleterious: not only will the increased dynamic range of the sensitivity function lead to a deterioration of the condition number of matrices to be inverted, but any experimental system will become highly sensitive to miss-location of the applied optical sources. Moreover, a greater amount of total optical power can be delivered to a tissue experimentally if it is illuminated by a diffuse large aperture source.
In this work we gained sufficiently independent sets of internal data to permit the simultaneous recovery of µ a and µ s by the use of multiple optical sources and detectors, but this may be undesirable in practice. It is not at present evident the way in which independence in the interior data can best be achieved (the choice of multiple optical source and detector locations signifies the implicit assumption that the independence is fundamentally related to the spatial gradient of the internal fields). This point is worthy of further consideration, since if independence could be achieve by another means, for example, an appropriate set of structured illumination patterns, this would have significant experimental advantages.

Measurement types
In this work we presented the use of the first-harmonic fluence as our measurement type. This data-type arises naturally in a number of UOT detection mechanisms. In applying these reconstruction techniques to experimental data the forward model will need to be normalised in some way to match the experimental observations. Like in the case of DOT, these 'coupling-coefficients' must remain constant through the experiment, lest significant error be introduced to the reconstruction.
UOT offers a more robust measurement type, that of the modulation depth, which in our linear formulation would be found as the ratio of the first-harmonic fluence to the total fluence. This measurement-type has the attractive advantage of being self-normalised, such that changes in the coupling coefficients, providing they are consistent across the power spectrum, will not affect the measurement. We employed a modulation depth measurement type in a lag-domain model presented in a previous work [17], and found the sensitivity functions to demonstrate improved spatial localisation, and complete suppression of the optical source and detector locations. Similar results were found experimentally by Gunadi and Leung [49]. This was manifested in improved reconstructions which were insensitive to perturbation near the source location, even when a point-source was employed. This is of importance in various clinical applications where superficial changes in blood volume during, for example, functional response, or therapy, may otherwise come to dominate an image. For these reasons, it would be of significant value to extend the analysis of this work to the modulation depth measurement type.

Summary
In this work we have provided an overview of forward modelling techniques in UOT. We have demonstrated that a computationally efficient frequency-domain linearised diffusion model can be found by the formal linearisation of a diffusion approximation to a correlation transport equation for UOT. We have derived and elucidated the form of the correlation measurement density functions which arise for the first-harmonic measurement type. These sensitivity functions were employed in the derivation of the gradient of the an objective function. Employing this gradient in a non-linear optimisation technique permitted the simultaneous reconstruction of images of the optical absorption and scattering coefficients in both two-and threedimensions.