Open Access
21 August 2019 High-density functional diffuse optical tomography based on frequency-domain measurements improves image quality and spatial resolution
Author Affiliations +

Measurements of dynamic near-infrared (NIR) light attenuation across the human head together with model-based image reconstruction algorithms allow the recovery of three-dimensional spatial brain activation maps. Previous studies using high-density diffuse optical tomography (HD-DOT) systems have reported improved image quality over sparse arrays. These HD-DOT systems incorporated multidistance overlapping continuous wave measurements that only recover differential intensity attenuation. We investigate the potential improvement in reconstructed image quality due to the additional incorporation of phase shift measurements, which reflect the time-of-flight of the measured NIR light, within the tomographic reconstruction from high-density measurements. To evaluate image reconstruction with and without the additional phase information, we simulated point spread functions across a whole-scalp field of view in 24 subject-specific anatomical models using an experimentally derived noise model. The addition of phase information improves the image quality by reducing localization error by up to 59% and effective resolution by up to 21% as compared to using the intensity attenuation measurements alone. Furthermore, we demonstrate that the phase data enable images to be resolved at deeper brain regions where intensity data fail, which is further supported by utilizing experimental data from a single subject measurement during a retinotopic experiment.



Functional neuroimaging is a valuable medical tool that is employed in a wide area of medical applications, with several noninvasive or minimally invasive neuromonitoring techniques for examining functional brain activity available for clinical use. Functional near-infrared spectroscopy (fNIRS)-based optical imaging is a noninvasive, relatively inexpensive technology applied in numerous applications where neuroimaging is crucial, such as functional brain mapping,1 psychology studies,2,3 intensive care unit patient monitoring,4 mental disease monitoring,5 and early dementia diagnosis.6 Recent advances in system design have also paved the way to promise wireless wearable systems, allowing neuroimaging to be performed in real-life contexts where it is otherwise not possible.7

fNIRS techniques rely on injecting light into the tissue, then measuring the emerging light at nearby locations. Light propagation in biological tissue is dominated by scattering events in the near-infrared (NIR) window, where low absorption allows measurements on the surface of the head to sample the brain cortex.8 Differential measurements at two NIR wavelengths, known as NIR spectroscopy (NIRS), can recover oxy (HbO2) and deoxy (HbR) hemoglobin concentration changes that convey information about functional brain vascular events. When more wavelengths are used, it is possible to additionally recover concentration levels of water, lipids, and even measures of metabolism such as cytochrome c.9

Due to the strong scattering of NIR in tissue, utilizing sparse optical measurements for neuroimaging leads to low spatial specificity as compared to fMRI. Imaging quality and quantitative accuracy can be significantly improved by arranging the NIR optical source and detector elements in a dense grid to provide overlapping measurements that allow tomographic reconstructions, a procedure known as diffuse optical tomography (DOT). Using a high-density (HD) DOT configuration that provides access to overlapping measurements at multiple source–detector (SD) separation distances has been shown to improve image quality10 to a point comparable to fMRI.1,11 In addition, instruments with a large dynamic range can utilize distant measurement neighborhoods of the HD array, which improves reconstruction quality and depth sensitivity,12 as sampling depth increases with SD distance.13

When continuous wave (CW) NIR light is used, measurements provide access to only light intensity attenuation from the tissue. By contrast, amplitude-modulated NIR light, typically in the range of 100 to 200 MHz and known as frequency-domain (FD) NIRS, enables additional measurements of a phase shift of the NIR light.14 Measurements of the phase shift are related to the average photon path length and can be used to directly estimate rather than assume scattering parameters within the tissue, which enables more accurate optical parameter recovery.15 FD-NIRS has been used in breast optical tomography for characterization and monitoring of lesions,16 joint imaging for arthritis detection and treatment monitoring,17 infant brain development monitoring,18 and brain trauma assessment.19 In previous studies utilizing FD NIRS, the phase measurements have typically been disregarded when recovering hemoglobin concentrations.2022 This is largely due to the fact that the accurate absolute measurement of the phase is not trivial and requires complicated and often time-consuming instrument calibration procedures, as the phase can be affected by a number of factors such as fiber length, fiber bending, and attachment angle, each of which can present long-term drifts and requires recalibration. The frequency-domain multidistance approach has been shown to provide increased sensitivity to cerebral oxygenation using intensity and phase data;23,24 however, the multidistance approach has been extremely challenging to adopt with an HD-DOT tomographic setup due to the strong demands on the dynamic range of the high bandwidth detector instrumentation.

Alternatively, measurements of picosecond modulated light (time-resolved or TR) provide a direct measure of the temporal point spread function (PSF) of the detected photons to account for both tissue absorption and scattering. However, systems with enough channels to provide overlapping measurements, thus achieving DOT, have a low temporal resolution (75 s for a full imaging cycle14), though various mathematical models have been proposed to increase the temporal resolution.25,26 Although TR systems are considered the golden standard within fNIRS, typically only spectroscopic measures are used to provide bulk tissue optical measurements, and these systems have not been widely adopted for HD-DOT.

The CW, FD, and TR measurements are all sensitive to hemodynamics and are therefore indirect measures of variations in brain function. An alternative method, known as the event-related optical signal or fast optical signal, uses just the phase measurements of FD systems and has been reported to be sensitive to neuronal activity.27 It has been suggested that this technique does not suffer from the hemodynamic delay effect and may be sensitive to scatter-related changes within firing neurons. However, in practice, this technique requires extensive stimulus repetitions to allow adequate signal averaging to obtain required signal-to-noise, thus canceling the real-time monitoring potential and resulting in experimental sessions of extremely long duration (60  min for one stimulus28) that are impractical for many applications.

To investigate hemodynamic changes in a practical setting, measurements of differential phase changes are much easier to apply and have shown high correlation with the fMRI blood oxygen level dependent signal, suggesting that phase changes convey information about functional hemodynamic phenomena.29 This is an expected result as changes in the absorption of a medium will not only affect the total number of detected photons (the measured intensity) but also the average photon diffusion path-length (the measured phase). How a change in absorption within the volume leads to differential changes in both (logarithmic) intensity and (linear) phase can be understood when considering the distribution of the time of flight of the photons. Early arriving photons travel a shorter distance and stochastically sample more superficial tissues while the later arriving photons travel longer distances and therefore tend to both stochastically sample deeper tissues and have a higher probability of undergoing absorption events.30 Therefore, the phase measurement-sensitivity profile to changes in absorption is deeper than the sensitivity of the intensity attenuation measurement because early and late photons contribute equally to the mean time of flight while the total intensity is dominated by early arriving photons.31

To date, FD measurements have not been used in the context of tomographic optical functional brain imaging and the effects of incorporating phase information together with intensity in HD-DOT has yet to be investigated. This work provides a comparison of imaging resolution and localization accuracy between CW and FD HD-DOT by analysis of simulated PSF over 24 subject-specific models including an experimentally measured noise model. The underlying theory is provided to demonstrate the origins of FD measurements. The benefits of FD HD-DOT as compared to CW HD-DOT are also demonstrated through an example of a focal activation map from a retinotopic experiment from a healthy subject.




Forward Problem

To model the light propagation in tissue, commonly known as the forward model, a numerical model based on the finite element method (FEM) has been employed to solve the diffusion approximation through the NIRFAST32 modeling and image reconstruction software package

Eq. (1)

where r defines a spatial location, ω is the modulation frequency, the diffusion coefficient κ=1/3(μa+μs) with μs being the reduced scattering coefficient, Φ is the photon fluence, μa is the absorption coefficient, cm is the speed of light in the medium, and q denotes a light source. High-resolution meshes (mean mesh element volume <1  mm3) have been used throughout this work, which improve the accuracy of the forward model at a significant increase in the computational cost.25 To overcome the computational cost, GPU parallelized linear solvers have been deployed that facilitate the efficient calculation of the forward problem in FD simulations; when using high-resolution meshes with 45,000 nodes, with 158 sources and 160 detectors, the forward model calculation time is less than 3 min.33


Inverse Problem

The sensitivity matrix, J, also known as the Jacobian or weight matrix is calculated using the Adjoint method34 that relates changes in boundary fluence measurements, y, with respect to small changes in underlying tissue optical parameters, x:

Eq. (2)


From Eq. (1), in the frequency domain, the measured data are complex numbers and can be decomposed in two parts: amplitude yI, and phase yΘ. Using the Rytov approximation, differential measurements with respect to a baseline y0 are described as a logarithmic ratio for amplitude and a difference in phase:

Eq. (3)


Therefore, the frequency-domain Jacobian describes changes in logarithmic intensity, lnyI, and linear differences in phase yθ, with respect to changes in the absorption, μa, and the diffusion, κ coefficients:

Eq. (4)


Assuming that only changes in the absorption coefficient are associated with vascular hemodynamic changes within the brain, this work will focus on modeling and reconstructing changes in absorption, μa, alone, reducing Eq. (4) to

Eq. (5)

and using measurements at two wavelengths of 690 and 830 nm, the dynamic changes of oxy (HbO2) and deoxy (HbR) hemoglobin concentration can be calculated using

Eq. (6)

where εc,λ is the extinction coefficient for chromophore c (HbO2/HbR) at wavelength λ. Conversely, it is possible to map the changes in oxy (HbO2) and deoxy (HbR) hemoglobin concentration directly to the measured dynamic data, by considering a spectral Jacobian:

Eq. (7)


This equation is used in conjunction with a noise model to produce realistic boundary data by perturbating hemoglobin concentration values.

To demonstrate the difference in the sensitivity distribution and magnitude of both log-amplitude and phase, as a function of source and detector distance, the Jacobian for a three-layered model, at a single wavelength of 830 nm, at a modulation frequency of 140 MHz is shown in Fig. 1. Here, the FD Jacobian has been calculated for four source and detector separations of 10, 20, 30, and 40 mm, for a three-layered two-dimensional model, with a skin/scalp layer of 5 mm, skull layer of 8 mm, and brain layer (combining gray and white matter) with optical properties shown in Table 1.

Fig. 1

Frequency-domain Jacobian for Log intensity (Joules/mm2/mM) and phase (radians/mM) for HbO2. Each row corresponds to a different source and detector distances of 10 to 40 mm. The dashed line represents the boundary of each layer, skin/scalp (top), bone (middle), and brain (bottom) tissue, each having different optical properties, as depicted in Table 1.


Table 1

Baseline optical properties in the head model.35

RegionOxyhemoglobin (mM)Deoxyhemoglobin (mM)Scatter powerScatter amplitude
Gray matter0.0560.0351.740.51
White matter0.0680.0271.310.82

Several features are evident from Fig. 1: (1) the Jacobian for the phase is showing sensitivity weighting at deeper regions as compared to log intensity, highlighting the potential of deeper tissue sampling, (2) the magnitude of change for log intensity measurements is larger than that compared to the phase but, (3) the distribution of the sensitivity for the phase appears more homogeneous across all tissues as compared to log intensity reducing the “hypersensitivity” to superficial signal contamination. These same features can also be demonstrated for HbR and is also true for a wavelength of 690 nm.

Functional DOT (fDOT) is concerned with three-dimensional spatial mapping of dynamic hemodynamic changes that occur during brain activations, using temporal data as measured from the scalp. Assuming that these vascular dynamic changes are small, the inverse model as used for parameter recovery can be considered as a linear problem. A Moore–Penrose pseudoinverse with Tikhonov regularization can be used to find an approximation of the inverse of the Jacobian in Eq. (2), J#, to perform a single step linear recovery of absorption:

Eq. (8)


Eq. (9)

where the regularization term is defined as α=a*max[diag(JJT)]. The weight of this regularization parameter, a, provides a tuning between high spatial frequency noise and smoothness in the image. Due to this association with image resolution,36 when running simulations without noise, fine-tuning of this parameter can improve image resolution beyond an experimentally and physiologically meaningful range. In practice, this parameter should be defined with respect to the noise of the measurements,10 therefore all simulations in this work use the same regularization weight of a=0.01, following previously reported methods.1

For the FD Jacobian, Eq. (5), the regularization is typically applied separately for each of the kernels of log intensity and phase measurements.32 In this work, the absorption coefficient is recovered for each wavelength then spectrally unmixed to recover HbO2 and HbR as described in Eq. (6). In previous work, spatially constrained reconstructions37 or spatially variant regularization schemes20,38 have been proposed that facilitate reconstructions of perturbations at deeper locations. To allow a direct evaluation of the potential increase in imaging performance afforded by the additional phase-based measurements, no spatially variant regularization methods have been considered in this work.


Subject-Specific Models

The subject-specific models are based on 14 female and 10 male individual subjects with a mean age of 26 (±4), with T1-weighted MPRAGE [echo time = 3.13 ms, repetition time (TR) = 2400 ms, flip angle=8  deg, 1×1×1  mm3 isotropic voxels] scans acquired for each subject. All subjects passed MR screening to ensure their safe participation. Informed consent was obtained and the research was approved by the Human Research Protection Office at Washington University School of Medicine.

High-resolution subject-specific FEM, consisting of five tissue layers were created, with mean element volume of 1  mm3, resulting in 450,000 nodes for each mesh using NIRFAST.39 A sparse-sensitivity matrix, J, (sensitivity threshold set to 0.1%40), was then created for each of the 24 subject-specific models with realistic optical properties for a 690- and 830-nm wavelengths, shown in Table 1. The reduced scattering coefficient was derived from the scatter amplitude and power, using μs=saλsp, where sa is the scatter amplitude, sp is the scatter power, and λ is the wavelength in micrometers.

In order to generate the boundary data, the nodes on the surface of the brain for each subject laying directly under the imaging array were selected to simulate point perturbations for HbO2 and HbR, Fig. 2(a). The data were then used for the calculation of a PSF to allow an analysis of the quality of the image reconstruction within a reasonable field of view. For each subject, the selected region depth from the surface of the brain varies, as the brain surface changes along the sulci and gyri morphology. To better demonstrate the depth variation of the surface of the brain, the histograms of brain surface depth across all 24 subjects are shown in Fig. 2(b) whereas Fig. 2(c) shows the brain surface depth probability distribution for each subject-specific model. This further highlights the importance and the need for improving the depth resolution of HD-DOT and demonstrates the need for modeling individual surface activation points for the presented analysis. The brain surface depth varies across the subjects, Fig. 2(c), ranging from 5 to 30 mm and brain surface points presented in this study lie on the cortical surface with a median depth of 12.5 mm, Figs. 2(b) and 2(c).

Fig. 2

(a) The modeled HD-DOT system with 158 sources (red) and 166 detectors (cyan), the surface of the brain surface is shown in pink and blue. The pink region denotes cortex locations of point perturbations in the simulations. (b) Histogram of brain surface depth across 24 subjects. (c) Brain surface depth probability distribution for each subject-specific model. The red stars mark the maximum probability for each subject while the red line is the most probable depth across all subjects (12.5 mm).



Realistic System Noise Model

The modeled system array adopted in this work is an HD-DOT system, with SD separation distances ranging from 13 to 47 mm as described previously.5 The SD pairs are extended over the whole of the head, in order to allow coverage of a large area of the optically accessible cortical surface of the brain that has previously been used as a benchmark for fDOT simulation studies.41 The modeled HD-DOT array contains a total of 158 light source at two wavelengths of 690 and 830 nm, and 166 detectors, resulting in 3,500 associated SD pairs within 50 mm separation at each wavelength, Fig. 2(a). The SD pairs are grouped into sets of nearest neighbor (NN) measurements (i.e., NN1, NN2, NN3, NN4, for first through fourth NN, respectively) with mean distances of 13, 29, 39, and 47 mm. In this work, the labeling for reconstruction with respect to a set of NN measurements used is the maximum NN; for example, for NN4 includes NN1, NN2, NN3, and NN4 measurements.

To generate a realistic noise model to be used with simulations, data were recorded using a frequency-domain system console (ISS Imagent™, Champaign, Illinois) employing six detectors at distances of 18 to 58 mm from the source, placed over the optical cortex. The source was modulated at 140 MHz and data were sampled at 39.74 Hz. Data were recorded for 120 s while the subject sat still and quietly fixated on a blank screen. Data were then processed following methods that have been previously described.1 Briefly, the data were band-pass filtered to 0.01 to 0.1 Hz to remove slow drifts, systemic physiology (heartbeat and respiration), downsampled to 1 Hz, and block averaged in 10 blocks of 12 s each. The same processing was applied to the measured amplitude and phase signals. The noise was then estimated as the standard deviation of the log mean intensity signals; therefore, it is a ratio-metric difference for the intensity and a linear phase difference (in degrees) for the phase and is shown in Fig. 3. As the system utilized photomultiplier tubes for signal detection, each detector was at a different voltage bias to increase signal gain and throughout the experiment, all data were ensured to be above the expected noise floor.

Fig. 3

In-vivo measured noise and a two-term exponential fit of noise versus distance from the source for intensity (left) and phase (right) measurements, on 690 nm (red diamonds) and 830 nm (pink squares) wavelengths.


The noise model n was calculated as a function of SD distance, r, via a two-term exponential of the form n(r)=aebr+cedr and was found to provide the best fit with coefficients shown in Table 2. In a previous work, using the same HD array configuration and similar signal processing but acquiring measurements with a CW system with LEDs at 830 nm, the amplitude noise was reported to be 0.12% for NN1, 0.15% for NN2, 0.41% for NN3, and 1.42% for NN4.12 The empirically derived noise model generated in this work corresponds to a noise level estimate for each neighborhood as shown in Table 3.

Table 2

The estimated coefficients for the two-term exponential fit of the noise model.

Amplitude noise modelPhase noise model
Coefficients830 nm690 nm830 nm690 nm

Table 3

The estimated noise levels for each measurement neighborhood.

Amplitude noise (%)Phase noise (deg)
NN830 nm690 nm830 nm690 nm

To simulate realistic measured data, perturbations of HbO2 and HbR values were induced using Eq. (7), using perturbation values scaled to minimize the difference between simulated and measured data. The perturbation for HbO2 was scaled to be 3.8  μM while for HbR was 1.8  μM. Evaluations were carried out using noise-free simulations and simulations that included the empirically derived noise model. Stochastic realistic noise was linearly added to each measurement based on sampling from a Gaussian distribution with zero mean and a standard deviation that was a function of SD distance with parameters given in Table 3.


Evaluation Metrics

To evaluate the image quality of each PSF, for each modeled focal activation, metrics as previously described in similar work have been utilized:10 (1) full width at half maximum (FWHM) defined as the maximum separation between all pairs of points in the activation above half of the maximum recovery; (2) localization error defined as the distance between the centroid of the recovery and the known location of the perturbed node; and (3) the effective resolution defined as the diameter of a sphere centered at the known perturbed node that can enclose all recovered nodes.

It is likely that the recovered volume above half of the maximum will not create a single focal recovery, resulting in multiple scattered activations.42 This is more evident in the cases with added noise, as the reconstruction often has multiple recovered activations. Therefore, to automate the quantification process, such activations are isolated by a connected component algorithm, using the FEM elements as the connectivity index and considering only the nodes that are associated with each other through at least one connected element. The recovered activations are then selected based on the maximum accumulated recovery value for evaluation without the use of any prior regional or locational knowledge. Therefore, a recovery near the perturbation will not necessarily be selected for evaluation, unless it is the strongest integrated recovery in the volume. Any recovery with localization error greater than 8 mm is considered a product of noise and therefore not included in the statistical analysis. The Wilcoxon signed rank test was used to assess the statistical difference between metrics of CW and FD reconstructions.


In-Vivo Experiment

An FD NIRS device (IMAGENT™, ISS Inc., Illinois) was used to obtain the NIRS-based data on a healthy subject. The system consists of 32 sources, modulated at 140 MHz, and 30 detector fibers with each source coupled to laser diodes emitting at 690 nm and 830 nm with the detectors being compact and fast photomultiplier tubes (PMT). Subjects were recruited from the University of Birmingham community, and written informed consent was obtained. Subjects were seated facing an adjustable screen while the imaging pad was attached over the occipital cortex with hook and loop strapping. Visual stimuli consisted of rotating logarithmic checkerboards, reversing intensity at 10 Hz, on a 50% intensity gray background. Each rotation completed in 36 seconds and 10 rotations occurred in the run, resulting in 7 minutes total scanning session time, including 30 second baseline periods at the beginning and the end. The raw signals were processed according to previously proposed pipeline,1 with all of the analysis being performed using NeuroDot, an extendable and user-friendly environment for analysis of optical data from raw light measurements (

Scans were in line with a previously reported CW high-density DOT imaging system setup.43 Briefly, the imaging cap consisted of 24 source positions (with two NIR wavelengths—690 nm and 830 nm—at each position) and 28 detectors interleaved in a high-density array. First- (13 mm, 1NN), through fourth NN (48 mm, 4NN) optode pairs were utilized, giving a total of 348 total possible measurements, at a frame rate of 40  Hz. Prior to data analysis, measurements with high noise (standard deviation of the logmean of the data larger than 7.5%1) were disregarded, resulting in 300 measurements passing the quality control for the scan. For the CW data, this was based on FD amplitude measurements only, whereas FD data contained both amplitude and phase. Image recovery of HbO2 and HbR were performed as defined above, using NIRFAST for the calculation of the Jacobian matrix for both CW and FD data.




Simulations Results

In this work, the brain surface under the imaging array was selected in 24 subject specific models and point activations modeled as perturbations were induced for each node in the selected area, allowing the recovery of PSFs both with and without noise. To provide an overview of the spatial distribution characteristics of the recovered PSFs, consider three points on the primary motor cortex within Brodmann Area 4, along the precentral gyrus, as shown in Fig. 4(a). These three focal activations are modeled at depths of 13.3 mm, 17.2 mm, and 20.5 mm, Fig. 4(b). PSF of each individual modeled focal activations for both HbO2 and HbR have been reconstructed and as both demonstrated similar results, for conciseness, only those of HbO2 are shown as in Fig. 5.

Fig. 4

(a) Axial view of the brain, the red line denotes the coronal slice selected to view the primary motor cortex, Brodmann Area 4, along the precentral gyrus. (b) The five layers of the model: skin, skull, cerebral fluids, white matter, and gray matter. Selected cortex points in different depths: 13.3 mm (left/red), 17.2 mm (middle/green), and 20.5 mm (bottom/blue).


Fig. 5

Recovered PSFs with different NN measurements (NN2, NN3, NN4), for activations at different depths. The perturbation location is denoted with the green dot. The FWHM contour (when present) is noted in cyan. Recoveries with (a)–(c) CW without noise, (d)–(f) CW with noise, (g)–(i) FD without noise, and (j)–(l) FD with noise).


For the focal point at 13.3 mm depth, Fig. 5, top row, both CW and FD provide qualitatively similar reconstructions for all SD pairs, but the FD demonstrated a visually better reconstruction in each case particularly for shorter NN pairs (NN1 and NN2) measurements. Specifically, for the noise-added CW recoveries, the localization error is 6.5 mm for NN2, 2.3 mm for NN3, and 1 mm for NN4, whereas for the FD recoveries are 2.2 mm, 1.6 mm, and 0.85 mm respectively. Similarly, the effective resolution for CW is 22.9 mm for NN2, 21.4 mm for NN3, and 11 mm for NN4, whereas for the FD recoveries is 17.7 mm, 16.2 mm and 9.8 mm respectively, demonstrating that for larger NN, the image recovery of the focal activation is improving for both CW and FD, and FD is consistently outperforming CW.

For a perturbation at 17.2 mm depth, Fig. 5, middle row, the noise added CW data fails to reconstruct for NN2 and NN3, whereas the NN4 has a 5.2 mm localization error and 17.9 mm effective resolution. The reconstructions based on FD data are qualitatively more accurate and have a localization error of 6.2 mm for NN2, 3.8 mm for NN3, and 2.2 mm for NN4; and the effective resolution is 21.42 mm, 20.58 mm, and 17.9 mm respectively. This case demonstrates that for some depth range the CW reconstructions may fail if the dynamic range of the detectors cannot provide measurements with a good signal to noise ratio (SNR), and the FD recoveries can perform reasonably well even when using only NN2.

When reconstructing for simulated data without noise, an activation positioned 20.5 mm deep from the surface of the head may be recovered reasonably well using any number of measurement neighborhoods, using CW only or FD data (Fig. 5) though the FD based data provide reconstructions with better image quality for all NN cases. In the presence of noise, recovery based on the NN2 and NN3 completely fail in both CW and FD. The recovery for CW using NN3 appears to provide a recovery near the perturbated point; however, this is not the strongest recovery in the volume, and therefore is considered as noise artifact and as such is not selected for evaluation. However, when using NN4, the FD data recovery demonstrates a better recovery, achieving 6.1 mm localization error with 24.8 mm effective resolution. In the presence of noise, recoveries for such deep points often fail, which leads to large increases in the localization and resolution metrics.

To provide a comprehensive evaluation of the benefits of FD data over CW, images were reconstructed for a set of NN1-NN4 SD pairs, using both noise-free and noise added data. In the presence of noise, the recoveries for deep activation points fail, which leads to significant increases in the calculated error metrics. The most obvious change is seen in the localization error, as beyond certain depths the recovered focal activations are not successful. In accordance to existing literature,42 any recovery with localization error greater than 8 mm is considered a product of noise. To identify the percentage of recoveries within the imposed localization limit, a success rate metric is introduced (Fig. 6), where a 50% success rate limit was set to examine the depths at which different approaches perform adequately. Recoveries exceeding this limit are not included in any further results or statistical analysis.

Fig. 6

HbO2 success rate for 8 mm localization error limit. The dashed blue and red lines denote the depth where the recoveries reach 50% success rate for CW and FD, respectively.


Each recovered focal activation is assessed using metrics of the FWHM (Fig. 7), localization error (Fig. 8), and effective resolution (Fig. 9). These results, quantized according to depth from the surface of the head, are summarized in Tables 4 and 5 and compared in Table 6. The noise-free evaluation is also presented to provide a benchmark of limitations imposed due to the ill-posed inverse problem, the employed regularization, and to provide a lower-bound on metric quality given further advancements in instrumentation that will lead to hardware with lower noise characteristics.

Fig. 7

HbO2 localization error median for CW (blue asterisk line) and FD (red cross line); the colored bands represent the 25th to 75th percentile range. The percentile bands stop at the point where the 50% success rate is reached, while the median line continues until the 10% success rate is reached.


Fig. 8

HbO2 FWHM median for CW (blue asterisk line) and FD (red cross line); the colored bands represent the 25th to 75th percentile range. The percentile bands stop at the point where the 50% success rate is reached, while the median line continues until the 10% success rate is reached.


Fig. 9

HbO2 effective resolution median for CW (blue asterisk line) and FD (red cross line); the colored bands represent the 25th to 75th percentile range. The percentile bands stop at the point where the 50% success rate is reached, while the median line continues until the 10% success rate is reached.


Table 4

Image quality metrics medians and quartiles for recoveries with CW noise-added model.

Depth (mm)FWHM (mm)Localization error (mm)Effective resolution (mm)
5 to 1011.610.3+13.7119.9+12.410.69.5+
10 to 1513.311.7+15.712.911.6+14.512.511.2+
15 to 2014.612.6+
20 to 2514.512.7+

Table 5

Image quality metrics medians and quartiles for recoveries with FD noise-added model.

Depth (mm)FWHM (mm)Localization error (mm)Effective resolution (mm)
5 to 1010.89.8+
10 to 1513.111.8+14.612.311.1+13.711.310.1+
15 to 2015.113.7+16.814.913.4+16.613.211.5+
20 to 251512.9+1715.313.4+17.414.812.7+

Table 6

FD reconstruction relative difference. The bolded cells are not statistically significant (p>0.01).

FWHM (%)Localization error (%)Effective resolution (%)
5 to 10777.551.532.220.716.1119.5
10 to
15 to 2042.49.647.244.837.617.915.714
20 to 2541.33.523.530.543.37.99.113

For the noise-free cases, the success rate maintains near 100% until it steadily declines at a characteristic distance (Fig. 6). The FD approach consistently performs better than the CW case in deeper regions. For CW recoveries, the 50% success rate is crossed at 17.18 mm for NN2, 22.27 mm for NN3, and 24.97 mm for NN4; and the FD 50% success rate points are at 23.04, 27.53, and 29.38 mm, respectively. For the noise-added cases, the success rate drops under 50% at shallower depths. Specifically, for the CW case, the 50% success rate is crossed at 14.87 mm for NN2, 16.59 mm for NN3, and 19.71 mm for NN4, and the FD 50% success rate points are at 18.43, 19.75, and 23.04 mm, respectively.

The FD data significantly improve the localization throughout the imaging domain for noiseless simulations (Fig. 7). At each depth, the use of FD data demonstrates a statistically significant improvement (p<0.01) in localization accuracy. For the noise-added case, the localization error is larger for activations deep within the brain, and the improvement of the FD reconstructions is again statistically significant (p<0.01) at each depth interval. Specifically, when comparing a depth range of 10 to 15 mm, the localization error for FD presents a 58.9% improvement for NN2 recoveries, 47.3% for NN3 recoveries, and 21.1% for NN4 (Table 6). These improvements are statistically significant for all modeled focal activations (p<0.01) using the Wilcoxon-signed rank test.

The FWHM for the noise-free cases increases with depth, but the FD reconstructions do not always significantly improve the FWHM (e.g., in the NN2 case); however, employing the FD approach results in slightly improved performance and smaller variations within each depth interval (Fig. 8). For the noise-added case, the FWHM still presents slight improvements when employing FD recoveries, however, the improvement is not always statistically significant (Table 6), particularly for NN2 at activation depths of greater than 15 mm.

Finally, the effective resolution with respect to depth from the head surface is shown in Fig. 9. As expected, in the noise-free cases, the effective resolution increases with depth. As the effective resolution is a function of both FWHM and localization, in the noise-added cases, the improvements observed are larger than those of the FWHM but smaller than the localization improvements. However, the observed improvement for FD data is statistically significant for each depth interval. Specifically, when comparing for a depth range of 10 to 15 mm, the effective resolution for FD presents a 20.6% improvement for NN2 recoveries, 13.7% for NN3 recoveries, and 10.5% for NN4 (Table 6). These improvements are all statistically significant for all modeled focal activations (p<0.01) using the Wilcoxon-signed rank test.

The results are summarized in Tables 4 and 5 where the reconstruction metrics for recoveries with noise-added data have been binned and analyzed in 5-mm depth intervals and the median, upper, and lower quartile for each interval is reported. Finally, in Table 6, the percentage relative change for reconstruction using FD data as compared to CW is shown.


In-Vivo Results

A single snapshot of a recovered focal activation for a single subject is shown in Fig. 10. Here the cross-sectional sagittal and transverse planes of a single focal activation within the left visual cortex are shown as a change in recovered HbO2. The CW reconstruction is shown to recover an activation at a much more superficial depth, mainly appearing in the area of the superficial tissue; whereas, the FD-based reconstruction is able to recover the activation overlapping with the cortex of the brain. It is important to highlight that the same regularization parameter for the inverse model has been used, namely ensuring that the improvement in image recovery and accuracy is primarily due to the inclusion of phase data within the FD-based reconstruction.

Fig. 10

In-vivo visual cortex activation (50% threshold, ΔHbO2 [a.u.]) using (a, b) CW data and (c, d) FD data.




To date, most HD-DOT imaging studies have concentrated on the utilization of CW data. Although some studies have reported the use of FD NIR systems for spectroscopic mapping of brain functions, these have been mainly limited to utilization of phase and amplitude data for the separation of optical scattering components from absorption. We have recently demonstrated the information content of phase data from a retinotopic experiment showing that the phase data do contain rich information, potentially arriving from deeper regions within the brain and therefore minimizing the superficial tissue contamination.44 In this work, we have undertaken a detailed study of the potential improvements to image quality gained through incorporating measurements of the phase along with intensity in the reconstruction process.

Utilizing a simple layered slab model, the Jacobean (i.e., the sensitivity of both log amplitude and absolute phase to absorption changes) has been calculated and shown in Fig. 1. Several features are evident: (1) both log-amplitude and phase data show a decrease in signal due to a change to absorption; (2) while the sensitivity for both log-amplitude and phase data increases in depth of sampling with SD separation, the phase data demonstrates deeper sensitivity for a given SD separation; and (3) the hypersensitivity in superficial tissues observed with log amplitude data is less pronounced for phase data. This supports our earlier findings that phase data are able to detect changes deeper within the brain and are less prone to be contaminated to changes within the superficial tissue that may mainly be attributed to systemic cardiac and respiratory signal.

Utilizing a set of 24 subject specific models, Fig. 2 highlights that the surface of the cortex can be from 5 to 30 mm deep with respect to the surface of the head, further highlighting the need for the utilization of data-types that will allow deeper sampling of the brain tissue. Although this can, in theory, be achieved through the use of larger SD separations, in reality, this is not possible due to system design and low SNR at measurement distances above 45 mm.1,43 To demonstrate that FD data-types (log intensity and phase) can provide a better sampling of deeper tissue, the subject-specific models have been utilized to simulate a set of data corresponding to focal activations at all cortex surface nodes lying underneath the imaging array, highlighted as pink regions, Fig. 2(a). In order to ensure that these simulated data are realistic, the noise characteristics of a commercially available FD system, ISS Imagent have been measured, as a function of SD-NN measurements (Fig. 3) and were added to the simulated data (Table 2). To investigate the benefits of using differing NN (NN2-NN4) measurements for activations at different depths, three sets of representative focal activation recovery maps are shown in Fig. 5. In all cases as the number NN measurements increases, the depth recovery improves. In each case, the presence of noise is shown to introduce artifacts within the recovered maps. Importantly, the utilization of FD data is shown to improve all metrics of image quality relative to using CW. In almost all cases, the localization accuracy of utilizing FD data is shown to be much better than those of CW, even including an activation at 20.5 mm where the CW data failed to recover a meaningful map.

To provide an overview of these benefits, images of focal activations for all noise-added simulated cortical activations have been reconstructed and analyzed to provide group average metrics of localization accuracy (Fig. 7) full-width-half-max (Fig. 8) and effective resolution (Fig. 9) with data summarized in Tables 46. For all metrics in the noise-free case, the accuracy decreases as a function of activation depth, with FD data outperforming CW in all cases. To achieve meaningful statistical comparisons, a localization error limit of 8 mm has been imposed as the upper-bound for statistical analysis, excluding all recoveries that are solely products of noise. For all NN measurements, the FD data provide a significant improvement in almost all recovered activations for all reported metrics (summary in Table 6). For example, assuming a set of NN3 measurements, it is shown that as compared to CW, the FD data show up to 59% improvement in localizing accuracy and up to 21% improvement in effective resolution. In addition, the FD data using NN4 are observed to successfully recover activations as deep as 20 to 25 mm, beyond the range accessible to CW.

Finally, to demonstrate the improvements using experimental data, a single focal activation from a set of retinotopic data for a single subject is shown in Fig. 10. Utilizing the same reconstruction algorithms, with the same regularization parameters, it is demonstrated that the inclusion of phase data from an FD system is capable of recovering the focal activation within the visual cortex with much more realistic location accuracy.

It is important to highlight that utilizing a time-resolved system for NIRS, utilization of mean-time of photon travel has been previously reported,45 which is analogous to the phase for an FD system. However in this work, the utilization of both log intensity and phase for tomographic reconstructions have been highlighted to demonstrate benefits. Although the sensitivity of the log intensity is shown to be more weighted toward the superficial tissues, it is this variation of sensitivity between log intensity and phase data that is providing the additional accuracy as demonstrated.



This work has examined the effect and potential benefits of employing frequency-domain measurements for functional DOT, using a high-density measurement array. Simulated data using 24 subject-specific models utilizing experimentally measured noise have been employed to provide a set of metrics for characterization of recovered image quality, as well as an in-vivo retinotopy experiment from a single subject. The results reveal that, as expected, resolution in the depth is related to the SD distance and level of the noise present, which are in agreement with previous studies.6 In addition, using the phase information significantly improves image quality when employing the same number of NN measurements with or without noise. While the use of the phase significantly improves image quality throughout the volume, the improvement is larger in the deeper regions, as phase samples deeper than amplitude for the same measurement. As such, the FD measurements can be used to resolve activations in deeper regions, where the CW completely fails. Overall, this study suggests that the use of FD measurements in the context of HD-fDOT expands the potential for improvements in image quality beyond the current state of the art.


The authors have no relevant financial interests in the manuscript.


This work has been funded by the National Institutes of Health (NIH) Grant Nos. R01EB009233-2, RO1-CA132750, NIMH-R21-109775, and NIHMH-K01-103594 and in part by Marie Sklodowska-Curie Innovative Training Networks (ITN-ETN) programme, under Grant Agreement No. 675332 (BitMap).



A. T. Eggebrecht et al., “Mapping distributed brain function and networks with diffuse optical tomography,” Nat. Photonics, 8 (6), 448 –454 (2014). NPAHBY 1749-4885 Google Scholar


Y. Minagawa, M. Xu and S. Morimoto, “Toward interactive social neuroscience: neuroimaging real-world interactions in various populations,” Jpn. Psychol. Res., 60 (4), 196 –224 (2018). Google Scholar


A.-C. Ehlis et al., “Near-infrared spectroscopy as a new tool for neurofeedback training: applications in psychiatry and methodological considerations,” Jpn. Psychol. Res., 60 (4), 225 –241 (2018). Google Scholar


W. Weigl et al., “Application of optical methods in the monitoring of traumatic brain injury: a review,” J. Cereb. Blood Flow Metab., 36 (11), 1825 –1843 (2016). Google Scholar


V. Kumar et al., “Functional near infra-red spectroscopy (fNIRS) in schizophrenia: a review,” Asian J. Psychiatr., 27 18 –31 (2018). Google Scholar


S. Yomota et al., “Investigation of a task for early detection of cognitive decline in elderly people using functional near-infrared spectroscopy (FNIRs) for diagnosis of Alzheimer’s disease,” Alzheimer’s Dement., 13 (7), P1381 –P1382 (2017). Google Scholar


P. Pinti et al., “A review on the use of wearable functional near-infrared spectroscopy in naturalistic environments,” Jpn. Psychol. Res., 60 347 –373 (2018). Google Scholar


F. Jobsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science, 198 (4323), 1264 –1267 (1977). SCIEAS 0036-8075 Google Scholar


A. Gibson and H. Dehghani, “Diffuse optical imaging,” Philos. Trans. R. Soc. A, 367 3055 –3072 (2009). PTRMAD 1364-503X Google Scholar


B. R. White and J. P. Culver, “Quantitative evaluation of high-density diffuse optical tomography: in vivo resolution and mapping performance,” J. Biomed. Opt., 15 (2), 026006 (2010). JBOPFO 1083-3668 Google Scholar


M. D. Wheelock, A. T. Eggebrecht and J. P. Culver, “High-density diffuse optical tomography for imaging human brain function,” Rev. Sci. Instrum., 90 051101 (2019). RSINAK 0034-6748 Google Scholar


H. Dehghani et al., “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Appl. Opt., 48 (10), D137 –D143 (2009). APOPAI 0003-6935 Google Scholar


M. Firbank, E. Okada and D. Delpy, “A theoretical study of the signal contribution of regions of the adult head to near infrared spectroscopy studies of visual evoked responses,” Neuroimage, 8 (1), 69 –78 (1998). NEIMEF 1053-8119 Google Scholar


R. J. Cooper et al., “MONSTIR II: a 32-channel, multispectral, time-resolved optical tomography system for neonatal brain imaging,” Rev. Sci. Instrum., 85 (5), 053105 (2014). RSINAK 0034-6748 Google Scholar


M. A. O. Leary et al., “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett., 20 (5), 426 –428 (1995). OPLEDP 0146-9592 Google Scholar


S. Srinivasan et al., “Near-infrared characterization of breast tumors in vivo using spectrally-constrained reconstruction,” Technol. Cancer Res. Treat., 4 (5), 513 –526 (2005). Google Scholar


A. H. Hielscher et al., “Frequency-domain optical tomographic imaging of arthritic finger joints,” IEEE Trans. Med. Imaging, 30 (10), 1725 –1736 (2011). ITMID4 0278-0062 Google Scholar


H. Bortfeld et al., “Assessment of infant brain development with frequency-domain near-infrared spectroscopy,” Pediatr. Res., 61 (5), 546 –551 (2007). PEREBL 0031-3998 Google Scholar


A. Ghosh et al., “Normobaric hyperoxia does not change optical scattering or pathlength but does increase oxidised cytochrome C oxidase concentration in patients with brain injury,” Adv. Exp. Med. Biol., 765 67 –72 (2013). AEMBAP 0065-2598 Google Scholar


J. P. Culver et al., “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab., 23 (8), 911 –924 (2003). Google Scholar


S. D. Power, T. H. Falk and T. Chau, “Classification of prefrontal activity due to mental arithmetic and music imagery using hidden Markov models and frequency domain near-infrared spectroscopy,” J. Neural Eng., 7 (2), 026002 (2010). 1741-2560 Google Scholar


R. Zimmermann et al., “Detection of motor execution using a hybrid fNIRS-biosignal BCI: a feasibility study,” J. Neuroeng. Rehabil., 10 (1), 4 (2013). Google Scholar


F. Scholkmann et al., “Cerebral hemodynamic and oxygenation changes induced by inner and heard speech: a study combining functional near-infrared spectroscopy and capnography,” J. Biomed. Opt., 19 (1), 017002 (2014). JBOPFO 1083-3668 Google Scholar


J. Choi et al., “Noninvasive determination of the optical properties of adult brain: near-infrared spectroscopy approach,” J. Biomed. Opt., 9 (1), 221 –229 (2004). JBOPFO 1083-3668 Google Scholar


S. Powell et al., “Real-time dynamic image reconstruction in time-domain diffuse optical tomography,” in Biomed. Opt., OM4C.5 (2016). Google Scholar


L. A. Dempsey et al., “Whole-head functional brain imaging of neonates at cot-side using time-resolved diffuse optical tomography,” Proc. SPIE, 9538 953818 (2015). PSISDG 0277-786X Google Scholar


G. Gratton et al., “Shades of gray matter: noninvasive optical images of human brain responses during visual stimulation,” Psychophysiology, 32 (5), 505 –509 (1995). PSPHAF 0048-5772 Google Scholar


G. Gratton and M. Fabiani, “The event-related optical signal (EROS) in visual cortex: replicability, consistency, localization, and resolution,” Psychophysiology, 40 561 –571 (2003). PSPHAF 0048-5772 Google Scholar


V. Toronov et al., “Study of local cerebral hemodynamics by frequency-domain near-infrared spectroscopy and correlation with simultaneously acquired functional magnetic resonance imaging,” Opt. Express, 9 (8), 417 –427 (2001). OPEXFF 1094-4087 Google Scholar


F. Martelli et al., “There’s plenty of light at the bottom: statistics of photon penetration depth in random media,” Sci. Rep., 6 1 –14 (2016). SRCEC3 2045-2322 Google Scholar


T. Binzoni et al., “Depth sensitivity of frequency domain optical measurements in diffusive media,” Biomed. Opt. Express, 8 (6), 2990 –3004 (2017). BOEICL 2156-7085 Google Scholar


H. Dehghani et al., “Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng., 25 (6), 711 –732 (2009). CANMER 0748-8025 Google Scholar


M. Doulgerakis et al., “Toward real-time diffuse optical tomography: accelerating light propagation modeling employing parallel computing on GPU and CPU,” J. Biomed. Opt., 22 (12), 125001 (2017). JBOPFO 1083-3668 Google Scholar


S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part 2: finite-element-method calculations,” Appl. Opt., 34 (34), 8026 –8037 (1995). APOPAI 0003-6935 Google Scholar


A. T. Eggebrecht et al., “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” Neuroimage, 61 (4), 1120 –1128 (2012). NEIMEF 1053-8119 Google Scholar


J. P. Culver et al., “Optimization of optode arrangements for diffuse optical tomography: a singular-value analysis,” Opt. Lett., 26 (10), 701 –703 (2001). OPLEDP 0146-9592 Google Scholar


D. A. Boas and A. M. Dale, “Simulation study of magnetic resonance imaging—guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt., 44 (10), 1957 –1968 (2005). APOPAI 0003-6935 Google Scholar


B. W. Pogue et al., “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt., 38 (13), 2950 –2961 (1999). APOPAI 0003-6935 Google Scholar


M. Jermyn et al., “Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography,” J. Biomed. Opt., 18 (8), 086007 (2013). JBOPFO 1083-3668 Google Scholar


M. Doulgerakis et al., “Towards real-time functional human brain imaging with diffuse optical tomography,” Proc. SPIE, 10412 1041209 (2017). PSISDG 0277-786X Google Scholar


X. Wu et al., “Evaluation of rigid registration methods for whole head imaging in diffuse optical tomography,” Neurophotonics, 2 (3), 035002 (2015). Google Scholar


Y. Zhan et al., “Image quality analysis of high-density diffuse optical tomography incorporating a subject-specific head model,” Front. Neuroenergetics, 4 6 (2012). Google Scholar


B. W. Zeff et al., “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U. S. A., 104 (29), 12169 –12174 (2007). Google Scholar


M. Doulgerakis, A. T. Eggebrecht and H. Dehghani, “Information rich phase content of frequency domain functional near infrared spectroscopy,” Proc. SPIE, 10865 108650C (2019). PSISDG 0277-786X Google Scholar


M. Diop and K. S. Lawrence, “Improving the depth sensitivity of time-resolved measurements by extracting the distribution of times-of-flight,” Biomed. Opt. Express, 4 (3), 447 –459 (2013). BOEICL 2156-7085 Google Scholar


Matthaios Doulgerakis is received his BEng degree in electronic engineering from Alexander Technological Educational Institute, Thessaloniki, Greece, in 2013, and his MRes in medical imaging from Kingston University, London, United Kingdom. He is a PhD candidate developing real-time algorithms for functional diffuse optical tomography (DOT) at the University of Birmingham, United Kingdom. His research interests are inverse problems and optimization algorithms in medical imaging, particularly in diffuse optics.

Adam T. Eggebrecht received his PhD in physics in 2009 and is currently an assistant professor in the Mallinckrodt Institute of Radiology at Washington University School of Medicine. He research is focused on designing and applying diffuse optical tomography (DOT) imaging methods in the human brain to understand mechanisms of typical and atypical functional brain development. He directs development of the NeuroDOT software package for modeling and analyses of diffuse optical measurements of brain function. He is currently funded by an NIH career development award to optimize and apply DOT to the study of altered brain function in children with autism. He is also funded to further develop DOT for studies of brain development in neonates, infants, and toddlers, as well as patients with implanted deep brain stimulators.

Hamid Dehghani is a professor of medical imaging in the School of Computer Science at the University of Birmingham, United Kingdom, and an OSA fellow. He has published over 100 peer-reviewed journal papers in the area of image reconstruction and numerical modeling and has a long and established track record in the development of biophotonics based techniques with specific applications in in-vivo optical imaging.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Matthaios Doulgerakis, Adam T. Eggebrecht, and Hamid Dehghani "High-density functional diffuse optical tomography based on frequency-domain measurements improves image quality and spatial resolution," Neurophotonics 6(3), 035007 (21 August 2019).
Received: 9 May 2019; Accepted: 30 July 2019; Published: 21 August 2019 Logo
Cited by 56 scholarly publications.

Image quality

Phase measurement

Data modeling

Near infrared

Tissue optics


Back to Top