|Home | About | Journals | Submit | Contact Us | Français|
We present 3D linear reconstructions of time-domain (TD) diffuse optical imaging differential data. We first compute the sensitivity matrix at different delay gates within the diffusion approximation for a homogeneous semi-infinite medium. The matrix is then inverted using spatially varying regularization. The performances of the method and the influence of a number of parameters are evaluated with simulated data and compared to continuous-wave (CW) imaging. In addition to the expected depth resolution provided by TD, we show improved lateral resolution and localization. The method is then applied to reconstructing phantom data consisting of an absorbing inclusion located at different depths within a scattering medium.
In complement of traditional techniques such as functional Magnetic Resonance Imaging, Diffuse Optical Tomography (DOT) is emerging as a low-cost and portable method for non-invasive cerebral imaging [1–3]. Based on the measurement of diffuse near-infrared light attenuation through the scalp, skull and brain, DOT assesses local optical absorption due essentially to the concentrations in oxy- and deoxy-hemoglobin. Modeling of light propagation through the head and solving of the inverse problem enable imaging of total hemoglobin concentration and oxygenation in the brain.
Most DOT instruments used in neuroscience are continuous wave (CW) systems, but time-domain (TD) technology is a recent promising alternative with advantages compensating for its increased cost and difficulty of implementation. These advantages include absolute characterization of tissue optical properties [4,5] (both absorption coefficient μa and reduced scattering coefficient μs′), depth resolution with single source-detector separation [6,7], and better sensitivity to cortical activation [8,9]. TD systems introduce short pulses of light into the tissue (up to tens to hundreds of picoseconds) and measure the temporal point spread function (TPSF) of photons exiting after propagation through the tissue.
One application of TD system is the determination of the absolute optical properties μa and μs′ of a tissue. One to several source-detector pairs are placed on the surface of the studied organ. The measured TPSFs [10,11], or moments of these TPSFs [12–14] are fit non-linearly with a model of light propagation in a homogeneous medium [4,11], a two-layer medium [15,16], a more realistic tissue-type segmentation of the head , or a complete voxelized volume in the case of tomographic imaging [18,19].
However, for functional studies, differential imaging between a baseline state and an activated state is sufficient, and a linear relation between changes in absorption and changes in intensity is generally assumed valid for the small absorption changes associated with brain function. In the backprojection method widely used for CW data reconstruction [20–23], the change in intensity between each source and detector is translated into a local change in absorption. A two-dimensional image is obtained by interpolation of all local absorption changes. The method is based on a modified Beer-Lambert law where a differential path factor (DPF) accounts for larger propagation length than source-detector separation . Hiraoka et al. extended the concept to inhomogeneous media , by introducing partial DPFs describing the optical pathlengths in different tissue types, hence taking into account the contribution from different layers to the change in attenuation. However, no depth information is available with single distance CW data, and both cerebral and superficial activations are projected on a single imaging plane.
This limitation of NIRS imaging in terms of depth resolution can be overcome with TD data, where depth information is contained in the photons’ time of flights. Steinbrink et al. applied the concept of partial DPFs to the time domain by introducing time-dependant mean partial pathlengths (TMPP) . They used a model of fifteen 2 mm thick layers, and performed Monte Carlo simulations to calculate the time-dependant sensitivity to each layer. They thus obtained a sensitivity matrix which they inverted by singular value decomposition. They could distinguish between extra and intra-cerebral signals during brain activation, with a single source-detector pair measurement. Liebert et al. used a similar method but computed the sensitivity of three moments of the TPSF – integrated intensity, mean time of flight, variance – to each layer . They showed depth-resolved time-course of local perfusion after dye bolus injection on healthy subjects and stroke patients . We used a simplified 3-layer model – scalp, skull, and brain – and showed that we could experimentally distinguish between superficial systemic signals and cerebral activation signals during a motor stimulus on a single source-detector pair . In all these studies, depth resolution has been shown, but no 3D imaging was implemented, since a single source-detector pair was used. The technique can be extended to 2D imaging with depth resolution by simple interpolation of all source-detector pair’s measurements . However, this approach limits the lateral resolution to approximately the source-detector separation.
In this paper, we describe an actual 3D linear reconstruction for differential imaging, based on inversion of the forward sensitivity matrix calculated in different delay gates, in a similar way that is sometimes implemented in CW DOT imaging . We discuss the performance and limitations of the reconstruction technique with simulated data. As already demonstrated in previous papers, we show that TD offers the ability to localize the depth of absorption contrast, which is not achievable with single distance CW data. Furthermore, we demonstrate improved lateral resolution and localization for TD compared to CW. We then apply the reconstruction technique to phantom data obtained with our time-gated system .
In this section, we describe the formalism implemented for the image reconstructions.
For computational simplicity, the studied medium is modeled as a volume of 10 cm by 10 cm laterally, and 3 cm in depth, divided into nvox voxels of 2.5 × 2.5 × 2.5 mm3 (see Fig. 1(a)). The background optical properties were set to μs′ = 10 cm−1 and μa = 0.1 cm−1, unless otherwise stated. The probe set on the surface of the medium has a square geometry of 4×4 sources and 3×3 detectors (source-detector separation 2.5 cm) as shown in Figs. 1(a) and 1(b). Only nearest neighbor source-detector pairs are taken into account for the reconstructions.
We model light propagation in the medium with the analytical solution of the diffusion equation for a semi-infinite homogeneous medium, calculated with the image source technique , and extrapolated-boundary condition . Figure 1(c) shows the resulting TPSF obtained for one source-detector pair, all source-detector pairs yielding identical TPSF’s since the medium is assumed semi-infinite and homogeneous. Since our experimental TD device described in Ref. 30 is based on a time-gated detection of the TPSF, we use in our simulations gated detection. The measurements consist in the intensity at nGates delay gates, integrated over the gate width wGate. The gates considered in our simulations are presented in Fig. 1(c). To comply with our experimental parameters , each gate width was set to 300 ps and the delay between two gates to 500 ps.
Differential images are obtained using the normalized Born approximation for a change in absorption: ΔI/I0 = A Δμa, where the sensitivity matrix A linearly relates the changes in the absorption coefficient to the changes in intensity. Δμa is the vector of absorption changes at each voxel, of length nvox, and ΔI/I0 is the vector of changes in the normalized measured intensity, of length the number of measurements nMeas = nSD nGates where nSD is the number of source-detector pairs and nGates the number of delay gates for each pair. In the time domain, the terms of the sensitivity matrix A can be computed by convolution of the direct and adjoint Green’s functions :
where G(rS, rD, τ) is the time domain Green’s function (GF) solution of the diffusion equation at delay τ for a source at position rS and a detector at position rD, rj is the position of the jth voxel, and rS(i), rD(i) and τi are respectively the source position, detector position and delay of the ith measurement.
In practice, we first calculated the GFs solutions of the diffusion equation at each voxel of the medium in the frequency domain , for each optode (source or detector), at 101 frequencies (0 to 20 GHz, frequency step 200 MHz). The frequency-domain solutions were then Fourier-transformed to yield the GFs in the time domain with a 25 ps time step. To obtain the sensitivity matrix for each source-detector pair at a specific delay τ, the source forward GF at time τ′ and the detector adjoint GF at time τ-τ′ were multiplied, the result summed for τ′ varying between 0 and τ, and then integrated over the width wGate of the detection gate.
Note that alternative methods to compute the sensitivity matrix could be used, in particular analytical solutions in the time-domain for a perturbation, like described in Ref.  for a transmission geometry.
Figure 2 shows profiles of the obtained sensitivity matrix for the 8 delay gates depicted in Fig. 1(c) (delays between 0.1 and 3.6 ns with 500 ps between two gates; wGate = 300 ps). Figure 2(a) shows profiles of the sensitivity matrix along a source-detector vertical plane, and Fig. 2(b) shows (x,y) profiles of the matrix at depth z = 0.75 cm. A logarithmic grayscale was used to allow visualization of the large dynamic range of sensitivities of all 8 gates. As expected, sensitivities at longer delays go deeper inside the medium, and also probe a wider lateral region. From this additional information relative to traditional CW sensitivity profiles, we can expect both depth resolution and improved lateral resolution and localization.
From the measurements y = ΔI/I0 and the computed sensitivity matrix A, the reconstructed image is calculated by inversion of the sensitivity matrix: = pAinvy, where pAinv is the pseudo-inverse of matrix A, computed with the following regularization:
CW sensitivity matrix and data were simulated with the same model, by integration of TD data over all time steps. CW and TD reconstructions will be compared in Section 3 to assess the improvement offered by TD imaging.
In this section, we assume a point-like change in absorption, simulate the corresponding measurement vector, and reconstruct the 3D map of the absorption changes. This gives us the imaging point spread function of our reconstruction algorithm and enables us to assess the performance of the method and compare it to a CW reconstruction scheme.
Figure 3 shows examples of reconstructions for a point-like inclusion located at the lateral position 1, between source and detector, as defined in Fig. 1(a), and with a depth z varying between 0.5 cm and 3 cm by steps of 0.5 cm. Both CW and TD reconstructions are presented. The rendered volume shows the contour of 80% of the maximum in the reconstructed absorption change. The following parameters were used: 5 gates starting at delays 0.6, 1.1, 1.6, 2.1, and 2.6 ns with a width wGate = 300 ps (gates 2 to 6 in Fig. 1(c)); β = 20; α = 10−3; signal to noise ratio = 100 at the peak of the TPSF.
The following general observations can be made: the CW data always reconstruct the inclusion at the same depth. If a traditional Tikhonov regularization is used instead (i.e. if A is used in place of B in Eq.1), the CW reconstruction is pulled towards the surface (data not shown), where the sensitivity is maximum . The effect of the spatially varying regularization matrix L is to force the reconstruction deeper under the surface. However, the reconstructed depth is unchanged with different actual depths for the CW data.
On the contrary, the TD method reconstructs the inclusion deeper as its actual depth increases. For the true inclusion at a depth of 1.5 and 2 cm, the reconstructed depth is actually slightly over-estimated, which results from the effect of the regularization matrix L. As the true inclusion gets deeper, the reconstructed depth becomes under-estimated (see inclusion at true depth 3 cm).
We also observe that the reconstructed volume at 80% of the maximum contrast is smaller for TD than for CW reconstructions, showing improved lateral resolution for this location of the inclusion.
More quantitative and systematic assessment of the effect of different parameters and of the improvement of TD over CW will be investigated in the following paragraphs.
The reconstruction performances were assessed by a number of parameters. The location of the center of mass (COM) of the reconstructed inclusion was calculated by taking into account all voxels with a contrast above 80% of the maximum contrast in absorption: rCOM = (Σi, Vox≥80% Riri)/(Σi, Vox≥80% Ri), where ri and Ri are respectively the position and absorption contrast of the ith voxel. We define the localization error as the distance between the true inclusion and the COM, both in depth and laterally. We call lateral resolution the contrast-weighted sum of the lateral distances to the COM, over all voxels with a contrast above 80% of the maximum contrast: Res = (Σi, Vox≥80% Ri|ρi − ρCOM|)/(Σi, Vox≥80% Ri), where ρi and ρCOM are the lateral positions of the ith voxel and the COM respectively.
We studied the influence of the number of gates included in the reconstruction. Figure 4(a) shows the evolution of the reconstructed depth (depth of the COM) as a function of the true depth of an inclusion located at lateral position 1, for different number of gates included (starting from gate 1 on Fig. 1(c)). The reconstruction improves, more strikingly for deep inclusions, as more late gates are included, up to 6 gates, after which the reconstruction does not improve anymore as we include more noisy data. The first gate only brings minor improvement for superficial inclusion (data not shown), and does not contribute to the data reconstruction for deep inclusions. We tried other combinations (data not shown), and found that the best one with our parameters was 5 gates every 500 ps from 0.6 ns to 2.6 ns.
With these parameters, the inclusion can be reconstructed with a depth error under 15 % down to approximately 2.5 cm. The reconstructed depths of deeper inclusions continue to increase, but with larger errors, e.g. an inclusion at 3 cm is reconstructed at 2.3 cm (an underestimation of almost 25%).
By penalizing more the voxels with higher sensitivity, the L matrix enables us to reconstruct the inclusion better than a simple Tikhonov regularization . However this matrix has to be thresholded so that regions far from a source-detector pair do not inappropriately receive larger weighting. The effect of this threshold is presented in Fig. 4(b), where the depth of the reconstructed COM is plotted vs. the true depth of the inclusion, for CW data and TD data with the 5-gate combination discussed above, and for a β factor of 10, 20, 50 and 100. With CW data, the absorption is reconstructed at a constant depth whatever its true depth, and this reconstructed depth simply increases as β increases. For TD data, a smaller threshold (large β) enables better reconstructions for deep inclusions (z > 2.5 cm), but also leads to an overestimated reconstructed depth for medium-deep inclusions (z around 1.5 to 2.5 cm). Moreover, for inclusions very close to the surface (under 1 cm) located just underneath a source or detector as in position 3 from Fig. 1(b), a large β also leads to strong artifacts (data not shown). A β = 20 was used in the following simulations, enabling a good compromise for reconstruction of different depths between 1 and 3 cm.
The evaluation of the baseline optical properties of a medium is subject to uncertainty . Therefore, we tested the influence of an error in the background optical properties on the reconstruction. We do not present data for this study, but enunciate general results we observed. An overestimation of the scattering coefficient (μ′s,recon > μ′s,true) by 20% led to an inclusion being reconstructed closer to the surface by approximately 1 to 3 mm, and vice versa for an underestimation of the scattering coefficient. We did not observe an effect of an over- or under-estimated absorption coefficient (by up to 50%) on the reconstructed depth. For these results, the true depth of the inclusion was varied between 0.5 cm and 3 cm.
One major advantage of TD data is the depth information provided by the time of flight of photons. In Fig. 4(a), we showed the evolution of the reconstructed depth as a function of the true depth for one particular lateral position. In Fig. 5, we present the evolution of the reconstructed depth as a function of the lateral position of the inclusion, for a 1 cm and 2 cm deep inclusion. In both cases, the depth error is small (within 20% at 1 cm, and about 5% away from the edges of the medium at 2 cm).
For the inclusion located at 2cm below the surface, it also worth noting that the reconstructed depth varies very little with the lateral position of the inclusion. This means that the reconstruction method has no “blind zone”, and the depth reconstruction performance remains good for any position of the inclusion. Note that we do not compare the depth error with that obtained by CW measurements as no depth information is provided without overlapping or multi-distance CW measurement.
In this section, we study the performance of the reconstruction method in terms of lateral localization and lateral resolution, and compare it with a CW reconstruction method. Figures 6 and and77 show the evolution of, respectively, the lateral localization error and the lateral resolution, for CW (left) and TD (right) reconstructions, at two different depths of the inclusion (1 cm, top, and 2 cm, bottom) as a function of the inclusion lateral position.
For both depths, TD shows reduced error in the lateral localization and better lateral resolution. Importantly we also note that both the lateral resolution and the lateral localization for TD reconstructions are more uniform with the lateral position of the inclusion relative to the probe than with a CW reconstruction.
These improvements of TD over CW for the lateral localization and resolution can be explained intuitively based on the sensitivity profiles presented in Fig. 2. For later delay gates, the sensitivity profile of a given source-detector probes deeper into the medium giving us depth resolution, but also probes a larger region laterally providing additional information to improve lateral localization and resolution.
We have shown in a previous study  that CW and TD systems yield similar contrast to noise ratio (CNR) for typical depths of cerebral activation (1 to 2 cm) and source-detector separations (2 to 3 cm) used in functional brain imaging, with CW even giving slightly better CNR. The intuitive explanation of this result is the following: even though TD data yield better contrast to deep inclusion by selecting photons which have traveled deep inside the medium, these measurements are also impeded by a much higher noise due to the low level of light at late delays. However the CNR of an image has to be evaluated in regards to other metrics of the image, in particular its resolution.
In this section, we varied the regularization parameter α between 10−4 and 10, both for CW and TD reconstructions, and studied the evolution of CNR and lateral resolution. Figure 8 is a parametric curve of the image CNR as a function of the lateral resolution. The true inclusion is located at position 1, and at a depth z = 1.5 cm. We compute the CNR as the average contrast to noise ratio for all voxels of contrast above 80% of the maximum contrast. The image noise was obtained by propagation of the measurement noise by: σx2= pAinvσy2pAinvT, where σx2 is the image covariance matrix.
The plot illustrates the trade-off between CNR and resolution: as the regularization is increased, the CNR of the reconstructed inclusion increases, but is counter-balanced by a worsening of the lateral resolution. We observe that for an identical CNR of the image, the lateral resolution is improved by TD reconstructions compared to CW. Similarly, at identical lateral resolution, TD reconstruction enables a much higher CNR of the image.
Our TD system, based on a Ti: Sapphire pulsed laser and an intensified CCD camera acting as a parallel time-gated detector, has been described in previous publications [8,29]. At each detector position, we use a bundle of 7 fibers of different lengths by increment of 10 cm, enabling simultaneous detection in 7 delay gates by steps of 500 ps. For functional brain imaging, we developed a flexible probe consisting of a 4×4 array of sources and 3×3 array of detectors in a square geometry each, with a source-detector separation of 2.5 cm (same geometry as used for the above simulations and presented in Fig. 1(a)).
The system was tested on a liquid phantom containing a spherical absorbing inclusion. One half of the probe was set over a tank (19×19×9 cm3) filled with a solution of intralipid and ink (estimated optical properties: μa = 0.14 cm−1, μs′ = 10 cm−1), containing a hollow glass sphere of diameter 15 mm. The sphere was filled with the same intralipid solution with 70 times higher ink concentration. Images were acquired for 30 seconds with the sphere located at position 2 (as defined on Fig. 1(b)), followed by 30 seconds of acquisition with the sphere outside of the field of view of the probe. This procedure was repeated 5 times to increase signal to noise ratio. We repeated the experiment for three different depths of the inclusion: top of the sphere at a depth of 9.5 mm, 18.5 mm and 32.5 mm below the surface.
We calculated the contrast of the image as the percentage change in intensity between the “off” state when the sphere is outside of the field of view, and the “on” state when it is under the probe. We averaged the contrasts obtained from the 5 successive on/off states and used these data to reconstruct an image. We used a regularization parameter α = 10−2 for the first two depths of the sphere. However, for the deepest inclusion, the contrast to noise ratio was very low, and we had to change the regularization parameter to 10−1 to obtain a reasonable reconstruction. Figure 9 shows the reconstructions for all 3 depths of the inclusion, as well as the true position of the top of the absorbing sphere.
The superficial inclusion is reconstructed at a depth of 7.5 mm, the second at a depth of 12.5 mm, and the deepest at 22.5 mm. In all cases the reconstruction is closer to the surface than the true inclusion, even for medium-deep inclusions, which differs from our simulations. Note that our simulations used a point absorption at a particular depth, while the experimental setup consists of a 15 mm diameter inclusion. Moreover our hypothesis of small perturbation breaks down as the solution in the glass sphere is strongly absorbing, hence the linear assumption is probably not valid anymore. In addition, we had to use a larger regularization parameter α than in the simulations because of the increased noise level of the experimental data. The effect of a higher regularization parameter is also to pull the reconstruction towards the surface (simulation data not shown). However, the formalism developed here enables reconstruction of experimental differential data at three different depths down to 3 cm of the true inclusion.
We presented a 3D reconstruction technique of differential TD gated data, and showed with simulated data that the method allows both depth resolution and improved lateral resolution and lateral localization compared to traditional CW reconstruction. With our typical experimental parameters, we also showed that only a limited number of gates is useful for optimal image reconstruction. The technique was used to reconstructed dynamic phantom images, where a spherical absorbing inclusion was embedded at various depths within a scattering intralipid solution. The technique presented here for gated measurements would be easily adaptable to moments of TPSF, provided that careful treatment of the covariance matrix is used in the reconstruction regularization.
We thank Anand Kumar for useful discussions on the spatially varying regularization. This work was supported by the NIH grants P41-RR14075, R01-EB002482 and NIBIB-EB00790.