1. Introduction

Time-domain diffuse optical spectroscopy (TD-DOS) measures the distribution of times-of-flight (DTOF) of photons propagating through the sample from a source to a detector [1,2]. DTOFs are obtained statistically using time-correlated single photon counting (TCSPC) [1,3] and several short (picosecond) laser pulses detected using a fast single photon avalanche diode (SPAD) [2,4,5]. TD-DOS has been widely applied for biosensing applications in diffuse optical imaging (DOI) and diffuse reflectance spectroscopy (DRS) [4–7]. In most biomedical applications for in vivo tissue sensing, near-infrared (NIR) light is used to provide functional (e.g., hemodynamics) and structural (e.g., cell size/density) information about the underlying biological tissue by quantifying the medium’s absorption coefficient ${\mathrm{\mu}_a}$ and reduced scattering coefficient $\mathrm{\mu}_s^{\prime}$ [8,9]. When acquired in the time-domain, experimental measurements are better able to decouple ${\mathrm{\mu}_a}$ and $\mathrm{\mu}_s^\mathrm{^{\prime}}$ while also allowing for depth discrimination of absorption changes, compared to continuous wave measurements [2,10].

To determine a medium’s optical properties from a measured DTOF, a theoretical (forward) solution to the time dependent photon diffusion equation [4,11] is used iteratively as an inverse model to fit measurements for known source-detector separations (SDS) and tissue geometries [1,4,11]. A Green’s function approach [12,13] is used to develop an analytical solution – the medium’s temporal point spread function (TPSF) – to an idealized delta-function input, $\delta ({t - {t_0}} )$ with ${t_0} = 0$ ns. To account for experimental temporal shapes of sources and response of detectors, the forward model is computed with a convolution of the TPSF and a measured Instrument Response Function (IRF) [4,11]. The inverse problem is known to be significantly dependent on both the shape and temporal position of the IRF [14,15]. Therefore, for accurate reconstructions, it is critical to carefully measure both the temporal shape and absolute temporal position of the IRF relative to the measured DTOF, for each wavelength and detector channel [14,15].

Because experimental DTOFs are fit with a forward model comprised of the convolution of the TPSF and the IRF, it is important to know when the injected pulse enters the medium (i.e., ${t_0}$ in Fig. 1(a)). The difficulty in measuring ${t_0}$ is schematically illustrated in Fig. 1. We note that the IRF is impacted by both its temporal position ${t_0}$ and shape. Each collection geometry can also produce an IRF with a different shape due to the inclusion of a thin diffusor (Fig. 1(a) vs 1b) or from a larger distribution of possible photon paths (Fig. 1(c)). When measured with the configuration shown in Fig. 1(a), ${t_0}$ is taken to be the peak or barycenter of the measured IRF [15]. In practice, the configuration shown in Fig. 1(b) is used, where a strongly scattering and attenuation layer is placed between the source and detecting fibers which introduces some small delay $\mathrm{\Delta }{t_1}$, in the barycenter of the IRF, relative to Fig. 1(a) [15]. However, in reflectance spectroscopy, custom probes are epoxied at the sample end leading to immovable relative configurations of the source and detector fibers. Thus, acquisition of the IRF requires collection in reflection geometry (Fig. 1(c)) and from an appropriate reflecting material [4,15]. This increases the IRF barycenter position by $\mathrm{\Delta }{t_2}$ relative to Fig. 1(a). In order to fit a measured DTOF, knowledge of $\mathrm{\Delta }{t_2}$ is needed, and thus requires careful calibration for quantitative reconstructions [15,16].

 

Fig. 1. Left column: Experimental configurations for measuring the instrument response function (IRF) in different geometries. Right column: temporal relationships between the measured IRF and the distribution of times-of-flight (DTOF) for photons from an arbitrary turbid medium for each corresponding experimental configuration. (a) The ideal (best possible) geometry for measuring the IRF – by directly coupling the source and detection fiber. ${t_0}$ represents the launch time of the incident photon pulse into the medium. (b) A practically used configuration – by introducing attenuating and diffusing layers between the source and detector fiber. This introduces a time delay $\mathrm{\Delta }{t_1}$ relative to the configuration in Fig. 1(a). (c) Reflectance configuration – by directing the incident pulse onto a surface and reflecting the incidt pulse into detector. This introduces a larger time delay $\mathrm{\Delta }{t_2}$, with $\mathrm{\Delta }{t_2} > \mathrm{\Delta }{t_1}$.

Download Full Size | PDF

Several approaches have been explored to analyze collected DTOFs in terms of the optical properties as well as to quantitate functional changes and depth sensitives while limiting the contribution from the IRF [17–20]. In addition to fitting the entire DTOF using a theoretical model, methods for reducing the dimensionality of the collected DTOF have been developed such as exploiting time-dependent mean partial pathlengths [21,22], integrated photon counts in various time windows [17,23], moment analysis [17,18,20], and Fourier components of the DTOF [24]. Moment analysis has particularly been shown as a promising technique to limit the errors arising from IRF measurements, as those discussed above, while also being sensitive to identifying depth-dependent absorption changes in media [17,18].

Here, we examine three approaches to reconstruct measured DTOFs without requiring direct experimental measurements of $\mathrm{\Delta }t$. 1) A constrained Monte-Carlo Diffusion Theory (MC-DT) approach that operates by shifting the DTOF and the convolved TPSF and IRF to peak at $t = 0$ while strongly restricting $\mathrm{\mu}_s^\mathrm{^{\prime}}$ values in the inverse model. The imposed constraint on $\mathrm{\mu}_s^\mathrm{^{\prime}}$ was achieved using a recently developed technique [25] that translates the measured differences of DTOF peak times at two different SDS from Monte-Carlo lookup tables (MCLUT) into $\mathrm{\mu}_s^\mathrm{^{\prime}}$. 2) A calibrated-DT approach that uses a reference sample to calculate $\mathrm{\Delta }t$ as described previously [4]. This calculated $\mathrm{\Delta }t$ is used to shift the IRF before inverse fitting for each target phantom. 3) A free-shift DT approach – here a third parameter ${t_s}$ is included to be freely fitted with $\mathrm{\mu}_s^\mathrm{^{\prime}}$ and ${\mu _a}$ in the inverse procedure. A time-variable is introduced into the DT-model to provide the TPSF for a delta-function input at ${t_s}$. Relative merits and drawbacks of each technique using performance metrics across a large set of tissue simulating phantoms are discussed.

2. Materials and methods

2.1 Hardware

The instrumentation used for experimental measurements is as described previously [25]. Briefly, a super-continuum laser (SC400, NKT Photonics, DK) with a pulse duration <100 ps was spectrally filtered using a band-pass filter (SuperK VARIA, NKT Photonics, DK) with a repetition rate of 40 MHz. Laser pulses were delivered to and collected from the sample by 400-$\mathrm{\mu}$m diameter optical fibers placed with center-center separation $\mathrm{\rho }$ and in contact with the sample surface. The reflectance was measured using a SPAD detector (PMD-050, MPD, IT) that was electronically coupled to a time-correlated single photonic counting (TCSPC) board (SPC-130, Becker & Hickl, DE). The IRF was obtained in reflectance geometry by reflecting the source from a mirror into the detecting fiber that was covered by a piece of paper, as described previously [25]. The full-width half maximum of the IRF for each detection channel and wavelength was measured to be less than 80 ps. A custom optical fiber probe was used that consisted of 4 colinear optical fibers (400 $\mathrm{\mu m}$, NA 0.22) to form 3 detection channels. Each channel had a SDS of 5, 10, and 15 mm (measured from one fiber, at the edge of the collinear array, set to be the source). Because our system used a single detector, the detecting fiber head was manually switched to select for any specific SDS.

2.2 Phantom tests

For the measurements, a cylindrical glass container (8-cm diameter and 8-cm height) was filled with an aqueous solution of 20% Intralipid (IL) (Sigma-Aldrich; MO, USA) and dried bovine hemoglobin (Hb) (H3760; Sigma-Aldrich; MO, USA). Optical characterization of $\mathrm{\mu}_s^\mathrm{^{\prime}}$ was taken from the average values of the intrinsic reduced scattering coefficient of 20% IL performed by multiple independent research groups [26,27]. The accuracy of this optical characterization was confirmed by a comparison between experimental data of pure IL solutions and Monte Carlo simulations [25]. The intrinsic absorption coefficient of Hb was calculated by measuring the transmittance at four concentrations using a spectrophotometer and then extracting the intrinsic absorption coefficient from the linear fit. The absorption coefficient of the medium was then calculated by the intrinsic absorption coefficient along with hemoglobin’s relative mass fraction to water and IL [26].

Two separate phantom sets with two different scattering levels were prepared by using 20 mL of 20% IL (set 1) and 40 mL of 20% IL (set 2) mixed with 750 mL of deionized water. The absorption coefficients in each set were independently varied by eight serial additions of ∼250 mg of dry bovine hemoglobin. For each laser wavelength used, a total of 18 different samples of optical properties, combining nine values of (per set) ${\mathrm{\mu}_a}$ with each set having one of two levels of $\mathrm{\mu}_s^\mathrm{^{\prime}}$. Four laser wavelengths were used (with center wavelengths of 650 nm, 700 nm, 750 nm and 800 nm ${\pm}\,5$ nm) to create 72 different phantom tissue models with known optical properties. Three repeated measurements were performed on each sample, at all available SDS ($\mathrm{\rho } = \; 5,\; 10,\; 15\; $ mm), and signals were acquired for 30 s per sample.

2.3 Data analysis

Collected DTOFs were analyzed using a known solution from diffusion theory to simulate the TPSF in a semi-infinite homogenous medium [4,28]:

As discussed in Fig. 1, the time ${t_0}$ represents the photon launching time, which is typically considered as the peak or barycenter of the measured IRF when the IRF is measured by directly coupling the source and detecting fiber. Typically, the use of a ${t_0}$ parameter in Eq. (1) is not needed because ${t_0}$ for IRF and the DTOF measurements are equal. Thus, convolution of the IRF with $R(t )$ puts $\tilde{R}(t )$ on the correct time scale, relative to the DTOF. In our experiments, the IRF was collected as shown in Fig. 1(c) (in a reflection geometry) since the fiber probe was configured for reflectance measurements. This introduced a shift $\mathrm{\Delta }t$ in the IRF time scale. Thus, without knowledge of those time-shifts, $\tilde{R}(t )$ could not be directly compared to measured DTOFs. We consider three different procedures to recover both $\mathrm{\mu}_s^\mathrm{^{\prime}}$ and ${\mu _a}$ using time-resolved reflectance obtained at one or more SDS, that do not depend on directly requiring $\mathrm{\Delta }t$ or ${t_0}$ for reconstructions and refer to them as (a) MC-DT (b) Calibrated-DT, and (c) Free-shift DT.

In the MC-DT approach, both $\tilde{R}(t )$ and the DTOF are shifted so that they peak at $t = 0$ ns. A constraint is then placed on $\mathrm{\mu}_s^{\prime}$ in the inverse model using a previously described approach to estimate $\mu _s^{\prime}$ from relative peak-time differences with MCLUTs [25]. The fitting bound for $\mu _s^{\prime}$ was calculated from the MCLUT using ${\mu _a} = 0.001$ and 0.5 cm^-1 as a lower and upper bound, respectively given a measured Δt_max at SDS of 5 and 10 mm. Such an approach does not require absolute time scales as it operates on the relative difference between peak arrival times at two SDS which is a primary requirement for this study.

The calibrated-DT approach indirectly calculates $\mathrm{\Delta }t$ by fitting $\tilde{R}(t )$ to a DTOF collected from a calibration reference phantom with known optical properties [4] for which $\mathrm{\Delta }t$ is considered a fit parameter that minimizes the least-square error between measurements and reconstructions. A $\mathrm{\Delta }t$ is estimated for each wavelength and SDS to shift the IRF time scale before fitting the DTOFs from target phantoms.

In the free-shift DT approach, Eq. (1) is modified by making the replacement $t\; \to t - {t_s}$ and ${t_s}$ is used as an additional parameter in the inverse model to optimize for along with $\mathrm{\mu}_s^\mathrm{^{\prime}}$ and ${\mu _a}$. This has the effect of introducing an additional parameter in the forward model that allows for $\tilde{R}({t - {t_s}} )$ time scale to be moved to appropriately match the DTOF.

Analysis for all three methods utilized the following procedure. Collected DTOFs and $\tilde{R}(t )$ were normalized by their maximum count rate, optical properties were then determined by fitting the DTOF with $\tilde{R}(t )$ utilizing a Levenberg-Marquardt procedure to minimize the least-square error [4,29]. The sample was considered to have a refractive index of 1.35 to match reported values for IL solutions [27] with an external index of refraction of 1.5 to match the glass optical fibers. The fitting range included all count rates higher than 60% of the peak value on rising edge of the DTOF and 0.1% on the tail. The range of optical properties in the fitting procedure was taken to be between $1 \le \mathrm{\mu}_\textrm{s}^\mathrm{^{\prime}} \le 60$ cm^-1 and $0.001 \le {\mathrm{\mu}_\textrm{a}} \le 0.5$ cm^-1. Twenty random start values in the described ranges were utilized in the inverse procedure, and the calculated ${\mu _a}$ and $\mu _s^\mathrm{^{\prime}}$ that showed the lowest least-square error were taken as the converged optical properties. A one-way ANOVA test was performed for the three analysis methods considering the recovered absorption and scattering coefficients separately at each SDS. Data was averaged across all experimental concentrations and wavelengths and a post-hoc Tukey honestly significant difference test was performed to indicate significant differences between analysis methods.

3. Results

Example fits for the three methods (a) MC-DT, (b) Calibrated DT, and (C) Free-shift DT are shown in Fig. 2 considering the same experimental DTOF and IRF with expected (true) optical properties of $\mu _s^{\prime} = 11.1\; \textrm{c}{\textrm{m}^{ - 1}}$ and ${\mu _a} = 0.19\; \textrm{c}{\textrm{m}^{ - 1}}$. In the MC-DT approach (Fig. 2(a)), the IRF, DTOF, and $\tilde{R}(t )$ are shifted to peak at $t = 0$. By constraining $\mu _s^{\prime}$, the MC-DT approach had the largest residuals in fitting (R² = 0.94) due to the poor fit at early times (before the peak). However, it often resulted in higher accuracy in the recovered optical properties. The calibrated-DT approach is shown in Fig. 2(b). $\mathrm{\Delta }t$ was calculated from a solid phantom [30,31] reference and used to shift the IRF time scale. The time scales of both the DTOF and IRF are translated so the IRF peaks at $t = 0$ as shown in Fig. 2(b). A main difference between the calibrated DT and free-shift DT is the time scales between Fig. 2(b) and 2(c). In Fig. 2(c), the time scale of the IRF and DTOF are the same as the TCSPC measurements. The difference between the IRF and DTOF peak times in Fig. 2(b) and 2(c) are due to not accounting for the $\mathrm{\Delta }t$ shift in Fig. 2(c).

 

Fig. 2. Example diffusion theory fits of an experimental DTOF using three different approaches: (a) MC-DT, (b) calibrated DT, and (c) Free-shift DT. The true optical properties are $\mu _s^{\prime} = 11.1\; \textrm{c}{\textrm{m}^{ - 1}}$ and ${\mu _a} = 0.19\; \textrm{c}{\textrm{m}^{ - 1}}$. The recovered optical properties and the sum of squares of the residuals, R², for each approach are listed in the figure.

Download Full Size | PDF

In Fig. 3, we show the recovered optical properties for the three methods across the entire set of phantoms measured, at a SDS of 15 mm and illumination at 700 nm. These trends (not shown) were similar for other wavelengths used. The top and bottom row in Fig. 3 show the 2.5% IL and 5% IL scattering levels, respectively. The left and right columns show the recovered ${\mu _a}$ and $\mu _s^{\prime}$, respectively, as functions of the true ${\mu _a}$ of the medium. In all three reconstructions, the recovered $\mu _s^{\prime}$ showed largest variability while all three approaches consistently tracked linear recovery of ${\mu _a}$. MC-DT showed the best consistency and accuracy for recovery of $\mu _s^{\prime}$, given the constraints placed on reconstructions. The calibrated-DT also consistently tracked (an unchanging) $\mu _s^{\prime}$ but showed a change in slope from true values in ${\mu _a}$. The free-shift DT method showed considerable crosstalk between ${t_s}$ and $\mu _s^{\prime}$ (not shown), but this minimally impacted recovery of ${\mu _a}$.

 

Fig. 3. Optical properties recovered at the two scattering levels for the 15 mm SDS channel at 700 nm for the three different analysis methods. Each row represents one scattering level for a single wavelength. The left column shows the measured ${\mu _a}$ against the true ${\mu _a}$ while the right column shows the measured ${\mu _s}^{\prime}$ against the true ${\mu _a}$ for the MC-DT, calibrated DT, and free-shift DT approaches. The MC-DT had the best linearity and recovery of absolute optical properties.

Download Full Size | PDF

In general, all of these approaches performed better at the higher scattering level (5% IL vs 2.5% IL) in accordance with diffusion theory serving as a better approximation in higher scattering media. Within each scattering level, similar trends and accuracy were observed at the four wavelengths used. In total, the absolute accuracy in recovered ${\mu _a}$ was better (< 4% points) at 750 nm than 650 nm. However, in the phantoms, the expected absorption values at 650 nm were roughly double those at 750 nm, thus impacting accuracy of DT based analysis. All approaches showed expected linear changes in absorption as shown in Fig. 3. The MC-DT and free-shift DT approaches increasingly underestimated ${\mu _a}$ for larger absorption values, whereas the calibrated DT approach systemically overestimated ${\mu _a}.$ Although the free-shift DT approach showed the largest errors in $\mu _s^{\prime}$, this did not appear to significantly affect the recovery of ${\mu _a}$. In most situations, the free-shift DT approach matched the accuracy of the MC-DT approach in recovering ${\mu _a}$ while being significantly more accurate than the calibrated DT approach.

Percent errors in the absorption coefficient were calculated as ${\delta _{{\mu _a}}} = 100\; \times |{\textrm{true} - \textrm{measured}} |\; /\; |{\textrm{true}} |$ and shown in Fig. 4 for the calibrated DT (sky blue), free-shift DT (orange), and MC-DT (navy) approaches. These values represent the mean ${\delta _{{\mu _a}}}$ across all phantoms and wavelengths used and error bars represent the standard deviation in ${\delta _{{\mu _a}}}$. Large errors were noted for short SDS ($\rho = 5$ mm), displaying the well-known limitation of DT at such distances. The Calibrated and Free-shift DT methods were not statistically different (p > 0.05) when $\rho = 5$ mm. However, all other groups were statistically different (p < 0.05) for each SDS when using a post-hoc Tukey HSD test. Errors were at nearly 10% for $\rho = 10\; \textrm{and}\; 15\; \textrm{mm}$ with the MC-DT, which is comparable to ranges reported previously [4]. Free-shift DT was less accurate than MC-DT with ${\delta _{{\mu _a}}}$ of 15-20% at $\rho = 10\; \textrm{and}\; 15\; \textrm{mm}.$ The calibrated DT approach showed the largest mean errors at 40-50%$.$ Surprisingly, the calibrated DT approach did not show improvements in recovery of absorption when increasing SDS between 10 and 15 mm.

 

Fig. 4. Calculated percent error in the recovered absorption coefficient (${\delta _{{\mu _a}}})$ using calibrated DT (sky blue), free-shift DT (orange) and MC-DT (navy). ${\delta _{{\mu _a}}}\; $ is shown for each experimental source-detector separation where error bars represent the standard deviation when averaging across all wavelengths and scattering levels. Large errors are seen when $\rho = 5\; \textrm{mm}$, however ${\delta _{{\mu _a}}}$ is ∼10% when $\rho = 10\; \textrm{and}\; 15\; \textrm{mm}$ using MC-DT and ∼15-20% using free-shift DT. The MC-DT approach was able to provide more accurate estimates at all SDS and wavelengths.

Download Full Size | PDF

As noted previously, all approaches performed better at the higher scattering level (5% IL). Although ${\delta _{{\mu _a}}}$ for the two scattering groups from the MC-DT approach were not significantly different, the calibrated DT approach performed better at the higher scattering level with the free-shift approach performing similarly at each scattering level as shown in Table 1. For $\rho = 5\; \textrm{mm}$, accurate recovery of optical coefficients was not possible from any of the approaches. Though somewhat of an improved accuracy was seen in the MC-DT approach when using lower wavelengths (higher scattering): ${\delta _{{\mu _a}}}$ = 26.1 ± 5.5 at 650 nm compared to ${\delta _{{\mu _a}}}$ = 67.75 ± 6.6 at 800nm. At $\rho = 10\,\textrm{and}\,15\,\textrm{mm}$, there was no significant difference between these wavelengths.

Table 1. Percent error in recovered absorption coefficients (${{\boldsymbol \delta }_{{{\boldsymbol \mu }_{\boldsymbol a}}}}$) at two scattering levels. ${{\boldsymbol \delta }_{{{\boldsymbol \mu }_{\boldsymbol a}}}}$ in the MC-DT approach had a small dependence on scattering level, whereas calibrated DT performed significantly better at the higher scattering level (5% IL).

View Table

Figure 5 shows ${\delta _{{\mu _s}^{\prime}}}$ averaged across all phantoms for the four wavelengths and two scattering levels, at each SDS using calibrated DT (sky blue), free-shift DT (orange), and MC-DT (navy). Recovered ${\delta _{{\mu _s}^{\prime}}}$ was less than 8% across all wavelengths and SDS using MC-DT with slight improvements at longer SDS: ${\delta _{{\mu _s}^{\prime}}} = 5.6\; \pm 1.5$ for $\rho = 5\; \textrm{mm}$ and ${\delta _{{\mu _s}^{\prime}}} = 3.5\; \pm 0.8$ for $\rho = 15\; \textrm{mm}$. No differences were observed amongst wavelengths within the two scattering levels.

 

Fig. 5. Calculated percent error in the recovered reduced scattering coefficient (${\delta _{{\mu _s}^{\prime}}})$ using calibrated DT (sky blue), free-shift DT (orange) and MC-DT (navy). ${\delta _{{\mu _s}^{\prime}}}\,$ is shown for each experimental source-detector separation where error bars represent the standard deviation when averaging across all wavelengths and scattering levels. The MC-DT approach was able to produce accurate estimates (< 10%) for all SDS, while the free-shift DT was able to produce estimates < 20% when $\rho = 10$ and 15 mm. Calibrated DT was only able to provide accurate estimates (< 10%) at $\rho = 15$ mm.

Download Full Size | PDF

On the other hand, the recovered scattering coefficient using calibrated DT only produced accurate measurements (< 10%) of $\mu _s^{\prime}$ at $\rho = 15\; \textrm{mm}$. The free-shift DT approach provide reasonable estimates (< 20%) at both $\rho = 10\; \textrm{and}\; 15\; \textrm{mm}$. The calibrated DT and MC-DT approach were not statistically different (p > 0.02) for $\rho = 15\; \textrm{mm},\; $ however, each of the remaining groups at each SDS were statistically different (p < 0.001). Recovered $\mu _s^{\prime}$ values were not directly correlated to estimates of ${\mathrm{\mu}_a}$, as the approaches produced different estimates of $\mu _s^{\prime}$ for similar values of ${\mathrm{\mu}_a}$. This is particularly highlighted in Fig. 3 where the free-shift DT and MC-DT approach could produce values of $\mu _s^{\prime}$ 20% apart, but still yield values of ${\mathrm{\mu}_a}$ to within 3% of each other. Additionally, the increase in the accuracy of $\mathrm{\mu}_s^{\prime}$ using the calibrated DT approach at 15 mm did not translate to any increased accuracy in recovery of ${\mathrm{\mu}_a}$, as seen in Fig. 4 and 5.

4. Discussion

In this work, we compared three methods: 1) MC-DT, 2) Calibrated DT, and 3) Free-shift DT to calculate both the ${\mu _s}^{\prime}$ and ${\mu _a}$ from experimentally collected DTOFs. These techniques worked without knowledge of the photon launch time ${t_0}$. The metric for comparisons was gauged using percent errors of recovered optical properties from DTOFs measured in phantoms relative to their true values. All three approaches could use a system IRF that could be measured in a configuration different from that used to acquire the DTOFs. In other words, all three approaches account for an uncertainty of having a different launch time of the incident photon pulse ${t_0}$ from the IRF. The calibrated and free-shift DT approaches accounted for this uncertainty with a calibration phantom or by introducing a third time-shift ${t_s}$ parameter in the inverse problem, respectively. The MC-DT approach shifted the IRF and DTOF time scale to peak at $t = 0$.

In Fig. 6, we show the nature of the ill-posed inverse problem in recovering the optical coefficients using the reduced chi-squared ($\chi _R^2$) maps for each method. The $\chi _R^2$ maps for each method (MC-DT: Fig. 6(a); calibrated DT: Fig. 6(b); free-shift DT: Fig. 6(c)) were computed for an experimentally measured DTOF from a phantom with expected optical properties $\mu _s^{\prime} = 10.2\; \textrm{c}{\textrm{m}^{ - 1}}$ and ${\mu _a} = 0.17\; \textrm{c}{\textrm{m}^{ - 1}}$ at $\rho = 15\; \textrm{mm}$ using the three analysis methods over a range of optical coefficients. For the MC-DT and calibrated DT approach, $\chi _R^2$ maps are shown as functions ${\mu _a}$ and ${\mu _s}^{\prime}$. While, for the free-shift DT approach (Fig. 6(c)), we show $\chi _R^2$ as a function of all three fitting parameters ${\mu _a}$, ${\mu _s}^{\prime}$, and ${t_s}$.

 

Fig. 6. $\chi _R^2$ distributions when fitting $\tilde{R}(t )$ to a single DTOF plotted for the three methods: (a) MC-DT, (b) Calibrated DT, and (c) Free-shift DT. A unique minimum is not observed in the MC-DT approach while the calibrated DT approach produces a well-defined minimum. Although the free-shift approach produced a minimum, there is considerable crosstalk between $\mu _s^{\prime}$ and ${t_s}$.

Download Full Size | PDF

In both Fig. 6(a) and Fig. 6(c), we can see the impact of seeking to optimize $\chi _R^2$ in the inverse space, when ${t_0}$ is unknown. In the MC-DT approach this manifests with $\chi _R^2$ showing no well-defined minimum making convergence highly dependent on the initial starting values (or constraints). In the free-shift approach, although a global minimum is present, it shows the inherent crosstalk between ${t_s}$ and ${\mu _s}^{\prime}$. It also shows why the crosstalk did not significantly impact recovery of ${\mu _a}$ (i.e., the minimum contour mainly runs along the ${t_s}$ and ${\mu _s}^{\prime}$ plane). On the other hand, knowing ${t_0}$ as with the calibrated DT approach gives a well-defined global minimum (Fig. 6(b)). Figure 6(a) also demonstrates why a constraint on ${\mu _s}^{\prime}$ provides stability and therefore improved performance of the MC-DT approach.

Alternatively, as is done with the calibration-DT approach, the temporal position of the IRF can be determined indirectly using reference standard of known optical properties [4]. Although the approach shows clear advantage in the optimization of the inverse problem (Fig. 6(b)) it is limited by the accuracy of the calibrated time shift. In our analyses, the calibrated-DT approach systemically overestimated the absorption coefficient which can be explained directly as a consequence of an inaccurate time-shift calibration. These inaccuracies can arise due to differences in the optical properties of the reference phantom, or the method used to determine baseline optical properties for each phantom. In our case, we used commercially given optical properties for the reference solid phantom [31] whereas the Intralipid phantoms were determined from averaged literature values and spectrophotometer measurements. Better results were reported previously for the calibrated-DT based approach [4]. Imposing additional spectral constraints on scattering coefficients when fitting calibrated time-shifts could help improve such estimates.

Introducing a time-shift as a third (free) parameter in the fitting process offers the simplest way to account for uncertainties in ${t_0}$. However, this comes at the cost of weakening the fit (Fig. 6(c)), increasing computational time and potentially having significant crosstalk between recovered ${\mu _s}^{\prime}$ and recovered time shifts [29,32]. In this study, we did not observe notable differences in convergence time across the three methods (all three methods converged < 100 ms). For the free-shift DT approach, we allowed ${t_s}$ to be shifted by a maximum of 1 ns in the modeled $\tilde{R}({t - {t_s}} )$. Strong crosstalk between ${t_s}$ and ${\mu _s}^{\prime}$ were observed that limited the accuracy of recovered ${\mu _s}^{\prime}$ using this approach. However, it would be possible to further constrain ${t_s}$ (e.g., for known experimental configurations) that could increase accuracy of recovered ${\mu _s}^{\prime}$.

As shown in Fig. 6(a), shifting the time scales of the IRF and DTOF to peak at 0 creates a highly ill-posed inverse problem. To overcome this, a constraint on ${\mu _s}^{\prime}$ was imposed using a previously reported technique [25]. This approach had the best accuracy in recovery of both optical coefficients. Although the performance here was good, some limitations of this technique are to be noted. The maxima of the DTOF has the highest photon count-rates, but the shot-noise contribution in detector-electronics is proportional to $\sqrt N $, where N is the count-rate [33]. Thus, determination of the peak time is prone to electronic noise and these could be exacerbated by having different count-rates for different SDS and/or transport coefficients. Further, timing resolutions caused by photon counting jitter, can also vary across detectors. For example, photomultiplier tubes provide timing resolutions larger than 100 ps [34] while superconducting nanowire single-photon detectors can be under 5 ps [35]. Here, SPAD detectors were used and signal acquisition used long integration times (of 30 s) which typically provided timing resolutions of tens of picoseconds [36]. As the peak-time differences measured were much larger than instrument resolution, the approach worked. However, choices of different detectors and timing electronics can significantly impact the MC-DT approach. A similar approach as the MC-DT used here, but based on diffusion theory and using photomultiplier tubes has been reported previously [37]. The MC-DT approach also assumes the DTOF peak-time is equally affected by the IRF at each detecting channel (for the 2 SDS used for the peak-time difference) and thus could impact its extension to multi-channel systems. Finally, the MC-DT approach does not scale well for larger ${\mu _s}^{\prime}$ as peak-time differences map a larger range of possible $\mu _s^{\prime}$ values which necessitates the use of short SDS channels that limit the sensitivity to only superficial layers [38].

The previous discussion focused on the impact of the temporal shift of the IRF relative to the real system response, and analysis methods to account for such uncertainty. In practice, there are additional sources of error that can impact each of the three approaches. We have used the well-established diffusion approximation to analyze our results, but DT can be severely limited at the SDS studied here and at modelling early arriving photons [29]. We note that reconstructions using DT are strongly influenced by the temporal fit range due to DT not having uniform validity across each photon arrival time. Additionally, the shape of the IRF in addition to the temporal position can greatly impact recovered optical properties. Non-ideal boundary conditions between the DT modelling and the experimental system could also impact results. Specifically, the Calibrated-DT approach had different boundary conditions when collecting measurements on the reference solid phantom (epoxy on resin) compared to the experimental water-based phantoms (epoxy slightly submerged in water). Finally, the influence of using shorter integration times and therefore increasing noise levels in data was not investigated. All of these potential sources of error can impact results.

5. Conclusion

We investigated three approaches to overcoming the difficulties in accurately measuring the Instrument Response Function (IRF) in time-resolved reflectance spectroscopy. Particular focus was given to overcoming the uncertainties in measuring the launch time ${t_0}$ of the incident photon pulse in reflectance geometries. We compared approaches that (a) operate by shifting the peaks of both the measured DTOFs and theoretical forward model after convolution with the IRF by overlapping them at $t\; $= 0, (b) calculate ${t_0}$ from a reference standard with known optical properties, and (c) introduce ${t_s}$ as a fitting parameter in the inverse problem.

We find that each approach has a set of unique advantages and shortcomings. When considering the three approaches and their respective inverse problem, having accurate knowledge of ${t_0}$ will provide the most well-posed reconstruction of optical properties (Fig. 6(b) vs 6a and 6c) leading to more accurate and reliable fitting procedures. As mentioned previously [29,32], the quantification of ${\mu _s}^{\prime}$ is more adversely impacted by uncertainties in ${t_0}$ than ${\mu _a}$. We show that recovery of absolute values of optical properties is still possible from time-domain reflectance using IRF measurements that are not exactly calibrated or that become uncalibrated in the presence of instrumental drift over the course of an experiment. Depending on the experimental setup and calibrations employed, accurate results can be achieved with all three methods to varying degrees.