Hierarchical Bayesian regularization of reconstructions for diffuse optical tomography using multiple priors

Abstract

Diffuse optical tomography (DOT) is a non-invasive brain imaging technique that uses low-levels of near-infrared light to measure optical absorption changes due to regional blood flow and blood oxygen saturation in the brain. By arranging light sources and detectors in a grid over the surface of the scalp, DOT studies attempt to spatially localize changes in oxy- and deoxy-hemoglobin in the brain that result from evoked brain activity during functional experiments. However, the reconstruction of accurate spatial images of hemoglobin changes from DOT data is an ill-posed linearized inverse problem, which requires model regularization to yield appropriate solutions. In this work, we describe and demonstrate the application of a parametric restricted maximum likelihood method (ReML) to incorporate multiple statistical priors into the recovery of optical images. This work is based on similar methods that have been applied to the inverse problem for magnetoencephalography (MEG). Herein, we discuss the adaptation of this model to DOT and demonstrate that this approach provides a means to objectively incorporate reconstruction constraints and demonstrate this approach through a series of simulated numerical examples.

Fig. 1

Diffuse optical imaging uses fiber optic based light sources and detectors to record changes in the optical absorption of underlying tissue. A grid of sensors is placed non-invasively on the head of a participant and used to measure changes in oxy- and deoxy-hemoglobin in the brain during task-evoked activation. The source-detector arrangement in the probe above is shown in Fig. 2.

Fig. 2

Simulation (A) and optical probe geometry (B) used in the construction of sample problems in this work. This probe was based on a tomography (over-lapping measurement) design described in [18] consisting of eight source positions (circles) and fifteen detector postions (squares).

Fig. 3

The optical inverse model was reparameterized in terms of wavelet coefficients. In the wavelet representation, the original image is described as a linear combination of low-pass, band-pass, and high-pass filter banks (left; for 1-dimensional case). The wavelet transform can be implemented in matrix form, which has the same filter structure (right) and will be used to apply a frequency bias to the superficial and deeper layers of the reconstruction.

Fig. 4

This figure shows a comparison the ReML and L-curve tuned Tikhonov regularized reconstructions for simulations at low noise (signal-to-noise ratio of 100:1). In the top row (row-A), the original target (A1), the EM-reconstructed image (A2) and the L-curve reconstructed image (A3) of oxy-hemoglobin ( + 1μM simulated) is shown. In the bottom row (row-B) the original and reconstructed images of deoxy-hemoglobin (−0.25μM simulated). Notably, the ReML and L-curve techniques are nearly identical for this trivial case of only a single regularization hyper-parameter (λ = Λ₁/Λ₂).

Fig. 5

This figure shows a comparison the EM and L-curve tuned Tikhonov regularized reconstructions for simulations at high noise (signal-to-noise ratio of 5:1). The definitions of the subplots are identical to Fig. 4.

Fig. 6

In this figure, we compare the value of the hyperparameter (λ) determined by the L-curve and ReML technique (REML λ = Λ₁/Λ₂) for simulations with a contrast-to-noise ranging from 1:10 to 100,000:1 (half decade intervals). The L-curve and ReML techniques agree closely over this range implying that the ReML technique performs as well as the L-curve method for the trivial example of a single covariance component.

Fig. 7

In this figure, a perturbation in oxy-hemoglobin only (row A) was simulated. No deoxy-hemogobin changes were simulated (row B). In the Tikhonov regularized inverse [Eq. (19)], which applies the same regularization factor to both the oxy- and deoxy-hemoglobin parameters, the L-curve technique (A3 and B3) gave a reliable estimate for oxy-hemoglobin, but this same level of regularization resulted in a very noisy deoxy-hemoglobin image. The ReML approach (A2 and B2) used separate hyperparameters to regularize the two hemoglobin species and resulted in close estimation of both images. Row B shows the deoxy-hemoglobin results.

Fig. 8

In this figure, measurements were simulated to have a signal-to-noise ratio of 2:1 at the 830nm wavelength and only 1:2 at the 690nm wavelength. The resulting image reconstructions obtained via the ReML regularization using separate covariance components for oxy- and deoxy-hemoglon (A3 and B3) was very noisy (as expected at this very low SNR). The noise in the images was reduced when a third covariance component modeling the negative-covariance between oxy- and deoxy-hemoglobin was also included (A2 and B2). Row A and B show the oxy- and deoxy-hemoglobin images respectively. Subplot A1 and B1 are the simulated (truth) images.

Fig. 9

In this figure, we compare reconstructions of the two-layered model. In rows A and B the superficial and deeper layers are shown. Only the oxy-hemoglobin results are shown. A perturbation was simulated only in deeper layer (B1). In A2 and B2, we show the reconstruction using a covariance component that spans both layers (akin to conventional Tikhonov regularization). Here, the same regularization is applied to both layers and the reconstructed image is heavily biased to the upper layer and underestimated. In A3 and B3, we show reconstructions using separate covariance components for the upper and lower level, which allows a total of four (2 layers x oxy- and deoxy-hemoglobin) hyperparameters to be estimated via the ReML method. This allows an empirically determined spatially distributed regularization of the model that results in correct placement of the reconstructed object in the bottom layer. This result is nearly identical to a cortically constrained reconstruction (B4) where the top layer is masked and only the bottom layer reconstructed.

Fig. 10

In this figure, we compare image reconstructions in the case of a two-layered model with non-zero noise structure ( + 1μM) in the superficial layer (A1). Row A shows the top layer and row B shows the deeper (“brain”) layer. Only oxy-hemoglobin results are shown. In A2 and B2, we show the results using the reconstruction using the ReML approach with covariance components for the two layers but without any frequency bias (e.g. σ₁ = σ₂ = 1 voxel). In A3 and B3, the reconstruction using a low-frequency bias in the top layer is shown (σ₁ = 2.2 voxels and σ₂ = 1 voxel). In B4, the reconstruction with a cortical-constraint is shown, which artificially pulls the superficial noise into the bottom layer and results in a grossly overestimated signal.

Fig. 11

In this figure, a two layer model with two perturbations placed either 6 voxels (40mm; row A) or 2 voxels (1.3 mm; row B) apart. Only the deeper layer and only the oxy-hemoglobin results are shown. In A2 and B2, we show the reconstructions using the ReML model without any specific region-of-interest priors. In A3 and B3, we show the reconstructions using a statistical region-of-interest prior. The arrows indicate the magnitude of the simulated values.

Fig. 12

Here, we show cross-sections of the reconstructions shown in Fig. 11.

Fig. 13

In this figure, we demonstrate the effects of using an incorrect region-of-interest prior. The simulated true image (A1; SNR = 10:1) had a single perturbation in the second layer. The top layer and deoxy-hemoglobin results are not shown. In A2, we show the reconstructed image without any region-of-interest priors. In B1, we show the reconstructed image using the correct region-of-interest as a prior (prior is outlined in black). Finally in B2, we show the reconstruction using an incorrectly placed region-of-interest prior (outlined in black). Using the incorrect prior produced nearly identical results to the case in which no prior was used demonstrating that the ReML method correctly assigned a near-zero weight to the incorrect prior.