Linear reconstruction of absorption perturbations in coherent ultrasound-modulated optical tomography

Abstract. Discrete form sensitivity functions are derived for a number of measurement types employed in autocorrelation-based ultrasound-modulated optical tomography. The Jacobian for a particular problem is constructed from the set of such sensitivity functions arising from a raster scan of a focused ultrasound field through a turbid medium. A linear reconstruction of an absorption perturbation is performed in a simulated difference data experiment, and the performance of the measurement types is compared under different degrees of added noise.


Introduction
The purpose of imaging and sensing techniques such as ultrasound-modulated optical tomography (UOT) is the production of images or spatially defined quantitative measurements of some parameter of interest. In the field of optical biomedical imaging, this is typically an image of the optical absorption coefficient within a diffuse medium at a particular wavelength. When multiple wavelengths are employed, the data may be combined via knowledge of the absorption spectra of various chromophores to form physiologically relevant images of, for example, oxygen saturation.

Sensitivity and Image Reconstruction in UOT
In UOT, the sensitivity of the system is defined approximately by the product of the optical sensitivity and acoustic amplitude in the system. The vastly different spatial resolutions of the acoustic and diffuse optical aspects of UOT encourage various approaches to the generation of images or spatially localized measurements, depending upon the nature of the optical configuration and the applied acoustic modulation.
The most basic approach is that of direct mapping; when the measurement inherently posses sufficient and controllable spatial selectivity, a measured datum may be assigned directly to the region of assumed sensitivity in an image. This approach is common in systems where the energy that is used to probe the system travels in straight lines without significant diffuse reflection from the interfaces or multiple scattering by a participating medium, such as in X-ray or ultrasound imaging. In a sensing application, the "image" is essentially a single region, but the implicit assumption of spatial resolution remains. Depending upon the application, the spatial resolution afforded by a focused or time gated ultrasound field is arguably sufficient to delineate regions of biological interest. If we neglect acoustic absorption and assume a constant ultrasound excitation amplitude probes each region in a given experiment, a UOT image produced via direct mapping records the scaled optical sensitivity in the medium. 1 In transmission mode configurations of a UOT experiment with a large optical étendue (such as those employing photorefractive, 2-4 photo-refractive polymer, 5 spectral-hole burning, 6 CCD, 7 digital holography, 8 or speckle contrast 9 based detection mechanisms) the optical sensitivity through a homogeneous medium transverse to the optical axis is relatively constant when contrasted with the exponentially varying sensitivity parallel to the axis of light collection. Since optical absorbers embedded in a turbid region reduce the optical sensitivity of the system in their vicinity, a transverse scan of the ultrasound probe can provide useful images of the embedded absorbers. Few demonstrations have been provided of more complicated absorption patterns in which shadowing of the optical sensitivity may distort the image, or of systems where an image is generated by a scan parallel to the optical axis in a high étendue system.
In low étendue systems that employ point-like source and/or detection mechanisms, the optical sensitivity varies significantly throughout the image in all optical configurations and scan directions irrespective of the homogeneity of the medium. In the context of UOT, this point was highlighted succinctly by the images produced by Lev and Sfez. 1,10 In this case, and in that of a high étendue system with the significant heterogeneity, the fundamental problem is that the optical sensitivity of the system is not uniform.
In the case of pure diffuse optical tomography (DOT), it is the very nonuniformity of the optical sensitivity which is used to generate images of the optical absorption. In DOT, it is recognized that the spatial resolution for a pair of point-like optical source and the detector is unsuitable for the direct image reconstruction, and instead a model-based inversion procedure is employed. In the most general sense, this approach employs a forward model to predict the spatial sensitivity of the measurement and this is used to drive the update of an image based upon the measured data. The DOT problem remains ill-posed, and a-priori information regarding the solution is typically imposed to constrain the reconstruction. Bratchenia et al. 11 previously demonstrated the recovery of the optical absorption coefficient in a turbid medium via a model-based inversion. In their work, the UOT signal was described in the frequency domain, and the acoustically driven modulation of the coherent input light was linearized in a diffusion-style model, similar to that presented by Allmaras and Bangerth. 12 A three-dimensional recovery of the optical absorption coefficient was performed using an iterative nonlinear optimization employing the Levenberg-Marquardt algorithm.
In this work, we take an alternative approach. We employ a time-domain model of the UOT signal which naturally provides many of the oft-employed measurement types in UOT, such as the modulation depth. Our forward model considers the nonlinearity of the acoustically driven decorrelation of the optical field, which allows us to examine in greater depth the form of the sensitivity functions which arise in the UOT experiment. In an approach common in DOT, we linearize our forward model in the optical parameters; this leads to a one-step difference data reconstruction which resolves an absorption perturbation from an assumed known background.
This work specifically considers the case of single fiber point-source and point-detector measurements of the field autocorrelation under monochromatic acoustic excitation.

Autocorrelation Measurement in UOT
In an autocorrelation-based UOT experiment, we measure the lag (τ) dependent autocorrelation of the flux exiting the medium at a particular location (r) on the boundary, yðr; τÞ, r ∈ δΩ in response to the insonification of the coherently illuminated medium. In practice, this value must be derived from the intensity autocorrelation function via the Seigert relation. 13 Neglecting an exponential decay in the correlation function due to Brownian motion, the data contain an oscillatory component at the acoustic frequency and its harmonics; this is depicted in Fig. 1. In this work, we consider ultrasound pressures of amplitudes small enough to ensure that only the fundamental frequency is significant; more details of the assumptions employed are presented later when we consider the forward model in Sec. 2.3.
We consider three measurement types: 1. yðrÞ DC ¼ yðr; 0Þ is equal to the continuous wave (CW or DC) intensity in DOT.
2. yðrÞ AC ¼ yðr; 0Þ − yðr; T∕2Þ is the maximum variation, or AC component, of the correlation flux over an acoustic cycle of length T.
3. mðrÞ ¼ 1 − yðr; T∕2Þ∕yðr; 0Þ is the modulation depth flux, which is defined as the quotient of the previous two quantities.

Linear Image Reconstruction
In this work, we will employ a linear (difference data) reconstruction technique. Although strictly limited to the recovery of qualitative information, 14 its simple formulation allows us to investigate many basic aspects of the UOTinversion, including its sensitivity to noise, and the effects of the a-priori information required to circumvent the ill-posedness of the underlying problem. Consider a turbid medium parameterized by its spatially varying absorption coefficient, xðrÞ ¼ μ a ðrÞ. The measured data yðr; τÞ are found by the application of the nonlinear operator F yðr; τÞ ¼ FðτÞ½xðrÞ: (1) The operator F incorporates the physics of the problem under the given parameterization, the boundary conditions, a given set of optical sources and detectors, and a particular ultrasound configuration.
In a linear reconstruction, we typically consider the change in a measurement Δy ¼ yðr; τÞ − y 0 ðr; τÞ brought about by a perturbation in the parameters ΔxðrÞ ¼ xðrÞ − x 0 ðrÞ. Expanding Eq. (1) in a Taylor series around x 0 and dropping the spatial and lag dependent notation momentarily and then discounting the higher-order terms and rearranging, we get ΔxðrÞ ≈ ðFðτÞ 0 ½x 0 Þ −1 Δyðr; τÞ: F 0 is the first-order Frećhet derivative of the forward model. 15 In the case of the modulation depth measurement Expanding this expression in a series around a baseline x 0 , discounting higher-order terms and rearranging Our task in this work is thus to develop an expression for F 0 . In this work, we will develop the discrete form of this operator which is the Jacobian matrix, J. 16 The inversion procedure of Eq. (3) (or its equivalent modulation depth formulation) is then given by

Forward Model
Sakadžić and Wang 17 previously derived a transport equation which describes the propagation of correlation in a turbid medium under insonification by a monochromatic acoustic field. The authors state a diffusion approximation to this transport equation which is valid under a number of assumptions regarding the nature of the acoustic excitation and optical properties. The principle limitation is that phase increments accrued at successive scattering sites are uncorrelated which requires that ðk a l tr Þ ≫ 1, where ðk a l tr Þ is the product of the acoustic wavenumber and the optical transport mean free path μ 0−1 s . Other limitations, such as the maximum permissible acoustic pressure amplitude, are introduced by the convenient approximations in the derivation. A more detailed analysis of the requisite approximations is presented in the original work, where the model is experimentally validated, demonstrating a good agreement in the region of validity of the standard optical diffusion approximation. The diffusion approximation reads where gÞ is the reduced scattering coefficient, g is the scattering anisotropy, Gðr; τÞ is the correlation fluence in the medium, q 0 is an isotropic source term, P 0 is the pressure amplitude of the acoustic field, k 0 is the optical wavenumber in vacuo, n 0 is the refractive index of the medium, k a is the acoustic wavenumber, ρ is the density of the medium, v a is the speed of sound in the medium, ω a is the angular acoustic frequency, η is the acousto-optic coefficient, S a is the amplitude of scatterer displacement relative to the acoustic field, and ϕ a is the relative phase of scatterer motion to the acoustic field.
The term hðrÞ, proportional to the square of the acoustic pressure amplitude, describes the acoustically driven decorrelation of the optical field. For a given medium, and within the region of validity of the model (where k a l tr ≫ 1), the term is also dependent upon the square of the acousto-optic coefficient (which takes a value of η ≈ 0.32 in water 18 ), the square of the optical wavenumber, and has an approximately inverse relationship to the acoustic frequency. The effect of this term is considered for a wide range of ðk a l tr Þ in Ref. 19.
This correlation diffusion approximation was originally derived under the conditions of isotropic scattering where g ¼ 0, though it was demonstrated to provide a fair approximation under the replacement of μ s with μ 0 s .

Sources, boundary conditions and detectors
A collimated source of coherent light is approximated by an isotropic point source located at a depth 1∕μ 0 s below the incident surface. We employ a modified Robin boundary condition Gðr; τÞ þ 2An · κ∇Gðr; τÞ ¼ 0; r ∈ δΩ; (9) wheren is the outward normal to the boundary at r, and A depends upon the refractive index mismatch across the boundary. 20 The outgoing correlation flux is given by yðr; τÞ ¼ −n · κ∇Gðr; τÞ; r ∈ δΩ: In the case of an indexed matched boundary where A ¼ 1, Eqs. (9) and (10) may be combined to give yðr; τÞ ¼ 1 2 Gðr; τÞ; r ∈ δΩ; A ¼ 1: We solve the correlation diffusion equation by the finite element method. [21][22][23] Equation (7) is multiplied by a test function which obeys the boundary conditions, and whose zeroth and first derivatives are integrable over the domain. The boundary conditions of Eq. (9) are incorporated by the subsequent integration by parts. The domain is subdivided into a mesh of nonoverlapping elements joined at N vertex nodes. On this mesh, we define a set of piecewise linear basis functions such that u i ðr j Þ ¼ δ ij for i; j ¼ 1; : : : ; N where r j located at the j'th vertex node. We subsequently approximate the solution Gðr; τÞ ≈ P N j u j ðrÞG j ðτÞ. Selecting the basis functions in the weak formulation to be the same as the mesh basis allows us to write the resulting linear system of equations: where A is the finite element system matrix. We express the parameters of the forward model, xðrÞ, the absorption-like decorrelation function hðrÞ, and the diffusion coefficient κ using the same basis functions such that, for example, μ a ðrÞ ¼ P k u μ;k ðrÞμ a;k . Consequently, and To make a measurement at a point r i Eq. (11) is implemented by a linear operator D T in which the appropriate elements of the row vector are initialized: We may now express a single measurement of Eq. (1) in a discrete form The subscript m indicates that this expression considers a single UOT measurement involving one source location r j;m , one detector location r i;m , and a single ultrasound pressure distribution pðr; τÞ m . For a single measurement, the linear inversion of Eq. (3) may thus be written

Correlation Measurement Density Functions and the Jacobian
In Eq. (16), the term ðD T m AðτÞ½x 0 −1 m q m Þ 0 represents the sensitivity of the measurement to a perturbation in the parameters of the forward model. In the context of DOT, Arridge generalized such sensitivity functions in the framework of photon-measurement density functions. 22,24 This work followed earlier derivations of the measurement sensitivity to perturbations in absorption by Schotland et al., 25 who termed this the photon hitting density, and others. [26][27][28] We now explore the correlation measurement density functions (CMDFs) which arise in correlation-based UOT. We assume that both the scattering coefficient μ s and the acoustic field distribution pðr; τÞ does not change between measurement of the baseline and the perturbed state in our linear reconstruction. Assuming that the source term is independent of absorption, we take the derivative of the forward model, expanding the derivative of the inverse system matrix 29 The system is parameterized by only the absorption coefficient of the medium, and we thus take the derivative of the system matrix with respect to each basis coefficient μ a;k We substitute this result back into Eq. (17) which provides an expression for the k'th element of the m'th CMDF. By definition, this also defines the k'th column of the m'th row of the Jacobian which we seek To transform this expression into a more useful form, we exploit the reciprocity of the correlation diffusion equation. Arridge expresses this reciprocity in the context of the diffusion equation 30 This expression states that the correlation flux across a point on the boundary of the medium due to an internal source of correlation fluence is equal to the correlation fluence measured at the same point in the medium due to an adjoint source of correlation flux across the same point on the boundary of the medium, when this source is scaled according to the appropriate measurement operator, here D. 24 In the discrete case, this reciprocity manifests itself in the symmetry of the finite element system matrix. To proceed, we apply the reciprocity expressed in Eq. (20) to the expression of the Jacobian in Eq. (19) where we have exploited the symmetry of the system matrix, its inverse and derivatives, and the properties of the transpose operation. Finally, we denote the solution to the forward problem GðτÞ f;m ¼ AðτÞ½x 0 −1 m q m , and the solution to the adjoint problem (in which the source term is the measurement operator) Here, we see that the Jacobian for a perturbation of absorption in the k'th node can be found by taking the inner product of the correlation fluence in the domain due to the (standard) source term, with the product of the basis function derivative term and the correlation fluence in the domain due to the adjoint source. We may thus compute the CMDF for a given ultrasound configuration by two computations of the forward model.

Modulation depth CMDFs and Jacobian
The CMDFs and Jacobian for the DC measurement type are given directly by Eq. (22). The AC sensitivity functions are straightforward and the modulation depth sensitivity functions are given by substituting the two expressions into Eq. (5), as follows:

Regularization
Although we expect an improvement in the localization of our UOT sensitivity functions over that of pure DOT, our present model still describes a diffusion process: thus the problem of reconstructing an internal parameter distribution from the data measured on the boundary remains ill-posed, and effectively underdetermined. The consequence of these considerations is that even if the sufficient measurements were available to form a square Jacobian, direct solution of Eq. (6) would fail to produce a stable solution, as noise in the measurements amplifies high-frequency oscillations in the resultant reconstruction.
In order to seek a useful and stable solution, it is conventional to apply some form of regularization to the problem. One commonly employed approach is to make the assumption that a "good" solution is one which balances the norm of the residual kðJðτÞ½x 0 ÞΔx g − Δyk with the norm of the solution kΔx g k, or its derivatives. This approach is that of the general form of Tikhonov regularization in which we replace the matrix Eq. (6) with a minimization problem For the underdetermined case, this minimization procedure is equivalent to solving 31 In the following section, we employ zeroth-order Tikhonov regularization such that L ¼ I. We select the regularization parameter λ according to the l-curve method, though the point of maximum curvature was chosen by the inspection when the algorithm failed to find the correct point.

Correlation Measurement Density Functions
In Sec. 2.4, we derived the CMDFs for the three measurement types introduced in Sec. 2.1. We will now visually inspect these sensitivity functions to gain insight into the sensitivity and spatial localization offered by these measurements.
A two-dimensional square domain of side 50 mm is located with its bottom left corner at ð0; −25Þ mm. The domain is assigned a uniform absorption coefficient μ a ¼ 0.01 mm −1 , reduced scattering coefficient μ 0 s ¼ 1 mm −1 , and isotropic scattering g ¼ 0. The refractive index of the domain is 1.4 and a matched boundary is employed. An ultrasonic field with a Gaussian profile of full width half maximum 2 mm and peak amplitude of 0.5 MPa is projected through the plane. An isotropic point source is placed at position (1,0) mm to approximate a collimated source incident perpendicular to the boundary at (0,0). A diffuse detector is located at (50,0) mm and integrates the outgoing flux according to a Gaussian profile with the full width half maximum of 0.1 mm. Figure 2(a) depicts the zero-lag (intensity) sensitivity function for a given measurement. The approximated collimated source is indicated with the inward arrow at the appropriate location, and the center of the diffuse detector with the outward arrow. This figure demonstrates the characteristic absorption sensitivity profile often seen in the literature of DOT. [24][25][26] The scale of the linear plot is truncated in the positive direction; regions close to the source and detector demonstrate extremely high sensitivity to absorption perturbations. This sensitivity function is insensitive to the acoustic field location: each row of the DC Jacobian which corresponds to a particular optical source and detector pair will be identical. Since in this work we consider only one such pair, a reconstruction using this information can only result in a scaled and regularized image of the same form as the sole unique sensitivity function depicted in Fig. 2(a). To produce a more detailed image would require that many more optical sources and detectors be employed.
The AC measurement absorption sensitivity is presented in Fig. 2(b). In this image, the acoustic field is scanned to a location of (20,10) mm, as indicated by the crossed circle. This figure demonstrates significant sensitivity in the region of the ultrasound, however, the region of sensitivity extends outward from the insonified region toward the source and detector locations. Figure 2(c) demonstrates the modulation depth sensitivity with the acoustic field scanned to the same location, (20,10) mm. This figure demonstrates that measurement type is insensitive to perturbations close to the source, detector, and on a smooth line connecting the two. Placing an optical absorber along this path of insensitivity will reduce the amount of correlated and uncorrelated light reaching the detector by an equal amount such that while the overall light level falls, the modulation depth remains constant. This argument is clearly supported by the analytical form of the sensitivity function given in Eq. (4).
The modulation depth sensitivity functions demonstrate high sensitivity in the region of the applied acoustic field. Moving away from the acoustic focus, across the path of zero sensitivity, we encounter a region of negative sensitivity where a perturbation in absorption will attenuate more correlated light than uncorrelated light. In Fig. 3, we demonstrate the modulation depth sensitivity as we scan the acoustic field through the medium, with the same color scale as Fig. 2(c). Note in particular that when the acoustic field is focused in the region of greatest optical sensitivity (as seen in the DC sensitivity of Fig. 2(a)), the sensitivity function appears to have the highest spatial localization, with two regions of negative sensitivity appearing in the profile. This may be problematic in a direct mapping approach as the absorption perturbations will be averaged over the entire region of sensitivity, but is automatically accounted for in the image reconstruction process.

Linear Reconstruction: Noise, Regularization, and Comparison With Direct Mapping
Having examined the form of the individual CMDFs which together form the Jacobian, we reconstruct an image under varying degrees of noise. The purpose of this exercise is twofold: first, the procedure will evaluate the functioning of the reconstruction process; second, the regularization required as a function of noise will be assessed. We perform the reconstruction on both the AC and modulation depth measurements. The reconstruction is performed on an homogeneous background with a regular set of absorption perturbations such that any aberrations in the reconstructed image will be evident. The background domain is the same as employed in Sec. 3.1. The perturbation domain is formed by adding an array of 16 circular perturbations of radius 2.5 mm with a modified absorption coefficient of 0.02 mm −1 . The perturbed measurements were generated using an independent mesh from the background measurements in order to prevent committing the inverse crime of generating measurements by the same model by which they are inverted. 32 To demonstrate an "ideal" reconstruction, the perturbed absorption coefficient was interpolated from its mesh to the baseline mesh upon which the Jacobian is computed, and the resultant image is displayed in Fig. 4.
The measurement set was generated by scanning an ultrasonic field with a Gaussian profile of full width half maximum 2 mm and peak amplitude of 0.5 MPa through the domain from ð5; −25Þ to (45,25) mm at 0.5 mm increments in the x and y. This resulted in a total of 17 × 21 ¼ 357 measurements for both the baseline and perturbed measurement. Figure 5(a) shows a noise-free reconstruction using the AC measurement type. This provides a reasonable qualitative reconstruction of the absorption perturbations of Fig. 4. This is especially true when it is considered that the square region enclosing each absorption perturbation is, in the ideal case, sampled by only 25 ultrasound scan locations. We note, however, that there are large aberrations near the source and detector where the AC CMDFs demonstrate considerable sensitivity. As shown in Fig. 5(b), adding 0.5% of proportional Gaussian noise degrades the reconstruction to the point where the perturbations are just recognizable. The increased regularization required to reconstruct this image has also suppressed the aberrations near to the source and detector.   Absorption perturbation Δμ a [mm −1 ] interpolated from the perturbation mesh to the baseline mesh upon which the Jacobian is computed. The baseline mesh does not include features of the perturbed geometry and hence the boundaries are seen to be slightly blurred.
The reconstructions of Fig. 6 demonstrates the modulation depth reconstruction in the case of measurements without noise, with 0.5%, and 1% Gaussian noise. In the noise-free reconstruction, an excellent approximation to the actual perturbation is generated both in a qualitative and quantitative sense. As suggested by the CMDFs displayed previously, the source and detector locations are completely suppressed in the resultant images. As the level of noise in the image is increased, the sensitivity off-axis is gradually suppressed by regularization. In the case of 1% noise, only those perturbations directly on on-axis are perceivable, though this measurement type appears to perform better than the equivalent AC reconstructions.
To demonstrate the improvement of the reconstructed images over that of a direct mapping approach, Fig. 7 illustrates the image formed by directly assigning the modulation depth difference measurement, ΔmðrÞ, to pixels located at the centers of the acoustic foci. In the absence of noise, the absorption perturbations along the axis of optical sensitivity ( Fig. 2(a)) can be distinctly located. The direct mapping approach demonstrates similar spatial resolution to the reconstructed images: this is to be expected since the sensitivity of the modulation depth measurement was previously demonstrated to be highly localized to the region of the acoustic focus. Without compensating for the spatially varying optical sensitivity, those absorbers offaxis are imperceptible. Moreover, the resulting unit-less image provides no quantitative information as would be required if we wished to attempt a recovery of a clinically relevant parameter. It is a common feature of many hybrid imaging modalities that the relationship of the data to the desired parameter is often unclear. We also note that in the case of image reconstruction, a variety of rigorous regularization strategies can be applied to stabilize the recovery of the underlying parameters under varying degrees of noise. In a direct mapping experiment, we may attempt to reduce the deleterious effects of the noise by some form of spatial filtering (for example, convolution with a smoothing kernel), but any such approach is somewhat ad-hoc.

Discussion and Conclusions
The sensitivity functions and associated reconstructed images presented in this work provide considerable insight into the potential of autocorrelation-based UOT. It has been shown that by moving the focus of an acoustic field in the medium, additional spatial resolution can be introduced into the UOT reconstruction: in the case of CW DOT, this would require that additional optical sources and detectors were employed, with associated cost and complexity.
Reconstructions based directly upon the AC measurement type are, in the absence of noise, capable of reproducing a fair representation of the underlying absorption perturbations. A disadvantage of this measurement type is the significant erroneous absorption indicated near the source and detector positions: this is a direct consequence of the extreme sensitivity of the measurement in these regions. If an absorption perturbation were placed near the source or detector, the influence of this measurement change would dominate the resulting image.
In a reconstruction employing both the AC and DC measurements, regularization would suppress the contributions from the