1. Introduction

Time-resolved reflectance spectroscopy (TRS) has the potential to directly provide bulk scattering and absorbing coefficients of a turbid medium without making assumptions about the medium’s composition or structure [1–4]. In TRS, picosecond laser pulses (fired at MHz rates) are injected into a turbid medium and the multiply scattered and attenuated diffusely reflected pulses are detected using a fast single photon counting photodiodes that are usually coupled to time-correlated single photon counting (TCSPC) boards to measure the photon distribution of times-of-flight (DTOF) [4]. These measurements are then quantified using theoretical (or numerical) approaches to extract the absorption (µ_a) and reduced scattering (µ_s’) coefficients of the medium [2–5]. In biological media, accurately and efficiently recovering the optical transport coefficients can help parametrize a variety of functional and structural properties of biomedical and clinical interest [6,7].

Uncoupling µ_a and µ_s’ using experimentally measured DTOF is known to be a difficult problem that needs careful measurements and calibrations [1–5,8–10]. The measured DTOF is a convolved response of the theoretical temporal point spread function (TPSF) and the instrument response function (IRF) of the experimental system [11]. The TPSF represents the (Green's function) response of the tissue medium to an incident Dirac delta pulse, while the IRF represents the finite temporal profile of the incident laser pulse measured as it propagates through the detection optics [11]. Extraction of the optical properties is typically done by iteratively convolving a theoretical TPSF with the experimentally measured IRF and fitting the convolved response via non-linear least squares to the measured DTOF [8,12,13]. Since the theoretical TPSF is computed from known optical properties, the process yields optical properties for the medium at convergence.

The above approach requires both having appropriate theoretical estimates to accurately model the TPSF, and accurate estimates of the system IRF. Since the IRF is a function of the laser source, fiber optics, and detection electronics used to deliver and collect signals, the IRF must be measured for all detection channels and wavelengths used experimentally [14,15]. Curve fitting of the DTOF using diffusion theory (DT) is usually the most common method for extracting optical properties of turbid media using TRS [16]. Although the process requires iterative reconvolution based reconstructions and are usually performed off-line, growth of computational resources could render these methods to operate in real-time [16]. However, applications of such models to analyze biological tissues shown large variance in extracted coefficients [1–3,15,17] particularly in scattering as reported recently [18]. Inverse methods using DT can also be computationally intensive and have been reported to show crosstalk between the derived optical coefficients [1,3,9,10]. Thus, a method such as the one we have presented could potentially facilitate improved fitting with DT.

Limiting crosstalk between recovered coefficients is critical for accurate estimation of optical properties, which in turn impact clinical and diagnostic utility of these methods. Different approaches have been described to constrain inverse fitting algorithms and increase quantitative accuracy of recovered optical parameters [19–21]. For example, measurements at multiple wavelengths were used to constrain spectral properties of the medium in order to improve the accuracy of recovered parameters [21]. Thus, methods that can extract (or even constrain) the optical coefficients of the medium can help DT based models in estimation of optical properties.

Here, we present a simple-to-use, heuristic technique that facilitates estimation of the µ_s’ of a medium directly from experimentally measured DTOFs. This method uses lookup tables constructed using MC simulations and translates measured peak-to-peak time differences (Δt) in DTOF obtained at two distinct SDS into µ_s’ of the medium. We describe implementation of the method and demonstrate its application to experimental data obtained from well-described tissue simulating phantoms [13,22,23].

2. Materials and methods

2.1 Instrumentation

Figure 1(a) shows the schematic of the experimental system, where the input pulse from a super-continuum laser (SC400, NKT Photonics, DK) was spectrally filtered via a band-pass filter (SuperK VARIA, NKT Photonics, Denmark) and coupled into one end of an optical (source) fiber (diameter=400 µm; NA=0.22; length=1 m). The distal end of the source fiber illuminated a medium of interest and formed the sensing head of a custom-fabricated optical probe (Gulf Photonics, FL). The sensing head had three detection channels (made from three optical fibers identical to the source fiber) which were epoxied at distances of 5, 10, and 15 mm from the source fiber – thus, our sensing probe had a fixed geometry. Reflectance measured from a selected channel (identified uniquely using the SDS) and was coupled into a single-photon avalanche diode (SPAD) detector (PMD-050, MPD, Italy), that was electronically coupled to a time-correlated single photon counting (TCSPC) board (SPC-130, Becker & Hickl, Germany). The laser repetition rate was 40 MHz and an electronic sync signal from the laser was used to trigger the TCSPC for signal acquisition.

 

Fig. 1. (a) Schematic of the time-resolved system used to obtain reflectance measurements. Dashed lines represent where the different detecting fibers connected to a detector (SPAD – single-photon avalanche diode; TCSPC – time-correlated single photon counter). (b) Time-resolved measurements from an Intralipid phantom at the three experimental SDS. (c) Magnified view of data in (b) to show peak-peak time differences observed between all SDS pairs. Δt_{5, 10} indicates the peak time difference between SDS of 5 and 10 mm (Δt_{10, 15}: for 10 and 15 mm SDS; Δt_{5, 15}: for 5 and 15 mm SDS). (d) Measured IRFs for the three SDS at 650 nm.

Download Full Size | PDF

Figure 1(b) shows representative time-resolved reflectance measurements at the three experimental SDS used. Figure 1(c) shows a magnified view of the DTOF peaks. The IRF was measured to ensure equality for each channel across all wavelengths by reflecting the source pulse into the detecting fiber using a mirror with a white piece of paper placed over the detecting fiber. Figure 1(d) shows the IRF measured at 650 nm for all three detector channels used. The average root mean square error (RMSE) calculated over three orders of magnitude for all IRFs measured across all channels and wavelengths used was 0.008. For comparison, repeated scans from a fixed target surface had RMSE of nearly 0.001, while DTOFs at differing absorption and scattering coefficients produce RMSE values larger than 0.1. Temporal stability of the peak position was also monitored over several hours of continuous acquisition and was shown to vary less than 0.7% per hour.

2.2 Monte Carlo simulations

Although DT can be used to model TPSF measurements at multiple SDS, it is known to be inaccurate for modeling reflectance at early times and small SDS [1,4,12,13,17]. Therefore, a previously developed time-resolved Monte Carlo (MC) model [24] was used to calculate the TPSF in a semi-infinite, homogenous geometry for 48 different tissue models. Each TPSF had distinct µ_s’ and µ_a values spanning 3–18 cm^-1 and 1.0 × 10⁻³⁻0.5 cm^-1, respectively. Tissue models were generated by permuting six different µ_s’ values (3, 6, 9, 12, 15, 18 cm^-1) with eight different µ_a values (1.0×10⁻³, 0.026, 0.050, 0.075, 0.18, 0.29, 0.39, 0.50 cm^-1). Absorption values were chosen on a log-scale to more accurately sample changes in absorption over several orders of magnitude.

Simulations were run to generate the TPSF (with photons being launched into and detected from the medium using optical fibers matching the experimental probe) at SDS of 5, 10 and 15 mm for all 48 tissue models. A temporal resolution of 2 ps was used for simulations at SDS of 5 mm, while the temporal resolution was 10 ps at SDS of 10 and 15 mm. Each MC run used 3×10⁸ photons, an anisotropy equal to 0.7 (using the Henyey-Greenstein phase function), an index of refraction equal to 1.35 (to match the IL phantom) [22,25,26], and fiber diameters and numerical apertures of 400 µm and 0.22, respectively, to match experimental measurements. Simulated TPSF data were resampled (using a cubic spline) to have a resolution of 0.1 ps. The peak-times calculated from the cubic spline using two different temporal resolutions in generating the TPSF were compared to ensure no errors were introduced by resampling.

Representative TPSFs from MC simulations are shown at two SDS in Fig. 2(a) (SDS = 5 mm) and Fig. 2(b) (SDS = 15 mm) for three different media (circles: µ_s’ = 6, µ_a = 0.18 cm^-1; squares: µ_s’ = 9, µ_a = 0.18 cm^-1; triangles: µ_s’ = 6, µ_a = 0.39 cm^-1). Vertical lines show t_max for each simulated reflectance (derived from the resampled data). These data show that in media with identical µ_s’ an increase in µ_a causes a decrease in t_max (by 7% at SDS of 5 mm and 13% for SDS of 15 mm for media shown by circles vs. triangles), while in media with identical µ_a an increase in µ_s’ causes an increase in t_max (by 11% at SDS of 5 mm and by 35% at SDS of 15 mm for media shown by circles vs. squares). These figures also show that at short SDS and low scattering, the absorption coefficient only weakly influenced t_max.

 

Fig. 2. Monte Carlo simulated TPSF for three different media at SDS of (a) 5 mm and (b) 15 mm. Simulated data are shown in symbols while resampled splines (see text) are the solid lines. Vertical lines show the t_max for each simulation. Interpolated peak time differences for SDS of 5 and 10 mm (Δt_{5, 10}) and for SDS of 5 and 15 mm (Δt_{5, 15}) are shown in (c) and (d), respectively. Note that the time-difference scales are different in Fig. 2(c) and Fig. 2(d) with longer SDS having longer Δt . Δt_{5, 15} computed from t_max shown in Fig. 2(a) and 2(b) are marked on 2(d) (symbol, shape and color identify corresponding media). The grey rectangle in Fig. 2(c) represents how the range of µ_s’ is confined fixed Δt (shown for Δt_{5, 10} = 50ps).

Download Full Size | PDF

The difference (Δt) between t_max at SDS pairs of 5-10 mm, 10-15 mm and 5-15 mm, were obtained for each of the 48 simulations. For each SDS pair, Δt values across the 48 simulated tissue models were linearly interpolated to create 2D lookup tables of Δt values for µ_a spanning 1.0 × 10⁻³ - 0.5 cm^-1 (in 1.0 × 10⁻³ cm^-1 increments) and µ_s’ spanning 3-18 cm^-1 (in 0.01 cm^-1 increments). Several interpolations methods were tested (e.g. linear, cubic, spline) to create lookup tables and showed lesser than 0.1% variation in generated tables. These computed Δt distributions are shown for two different SDS pairs – in Fig. 2(c) (SDS pair of 5-10 mm) and Fig. 2(d) (SDS pair of 5-15 mm). Symbols for Δt_5,15 values computed for data shown in Fig. 2(a) and 2(b) are marked on the Fig. 2(d) (shape and color of the marker identifies the optical properties). The horizonal grey line in Fig. 2(c) illustrates how a given Δt value (calculated using SDS of 5 and 10 mm) confines the range of µ_s’ (shown by the shaded rectangle). In other words, Δt_5,10 = 50 ps could represent media with µ_s’ ranging between 7.1–10.7 cm^-1 and µ_a varying from 1.0×10⁻³ and 0.5 cm^-1.

2.3 Lookup method

Data shown in Fig. 2(c) and 2(d) form the crux of the lookup method. In the forward direction, specific values for µ_s’ (x-axis) and µ_a (color) identifies a single point in the shaded region – the ordinate of this point is Δt. In the inverse sense, given a Δt value (using a given SDS pair) it may be associated with several pairs of µ_s’ and µ_a values (identified by the horizontal line within the shaded figure in Fig. 2(c)). Thus, translation of a measured Δt for a pair of SDS into µ_s’ would require knowledge of µ_a, or vice-versa. As is evident from these figures, for fixed Δt, the limits on µ_s’ are much more restrictive than the limits on µ_a. For the range of optical properties considered here, the required Δt for a SDS pair could be produced by media with µ_a that span the full range of simulated values (1.0 × 10⁻³⁻0.5 cm^-1) but had a much more confined range for possible µ_s’ values. In other words, the uncertainty in estimation of µ_s’ given an incorrect guess of µ_a would be much lower than the uncertainty in estimating µ_a given an incorrect guess of µ_s’. Thus, our lookup method translates the measured Δt (from a given SDS pair) directly into µ_s’ via the interpolated 2D MC lookup table for the corresponding SDS, but needs an input value for the µ_a of the medium to function.

2.4 Estimating the absorption coefficient

Two different approaches were used to estimate the µ_a required for using the lookup method. The first approach, referred to here as the ‘tail method’, uses a limit derived from DT where the slope of the natural logarithm of the DTOF at long time scales is assumed to be directly proportional to µ_a (Eq. (1)) [27]. $R({r,t} )$ shown in Eq. (1) represents the measured DTOF and for longer times, the value of the estimated absorption from Eq. (1) would improve. Using this technique to accurately estimate µ_a would require deconvolving the IRF from the measured DTOF and having the reflectance signal span many nanoseconds. In practice, a limited dynamic range in the detection system and the drawbacks in deconvolving the IRF from the DTOF reduces the accuracy of this approach in estimating µ_a.

However, as we expect Δt to only weakly depend on µ_a, the tail method provides a simple and direct way to extract (even if crudely) a value for use with the lookup method to estimate µ_s’. Further, due to a limited dynamic range of our instrumentation, the DTOF time scale was shifted such that the peak occurred at $t$=0, and Eq. (1) was applied to analyze the measured reflectance signal from its peak value till it fell to 0.1% of the peak value. Derivatives were numerically calculated between successive times, smoothed and translated to calculate µ_a for each measured DTOF using Eq. (1) [27]. Because calculated numerical values were noisy and showed variations for different SDS (in the same media), µ_a was estimated as the averaged tail-method estimate for every scan, at all three SDS used.

In the second approach, referred to here as the ‘transmittance method’, the intrinsic absorption coefficient of hemoglobin was calculated by measuring the collimated transmittance using a spectrophotometer at the varying concentrations of bovine hemoglobin used experimentally (described in Sec. 2.5). µ_a was then calculated from the percent transmittance using Beer’s Law for each concentration, and the intrinsic absorption coefficient (${\epsilon _a}$) was determined from the slope of the linear fit of the calculated µ_a at each concentration. Expected (true) absorption coefficients of the medium distinct from values obtained using Eq. (1) were then calculated using ${\epsilon _a}$ and the mass concentration of absorber used. Dissolving large amounts of solid hemoglobin (∼255 mg for each concentration change) as well as uncertainties in cuvette calibrations in transmittance measurements could impact the accuracy of the estimated µ_a by this method.

2.5 Phantom preparation and measurements

Liquid phantoms were prepared to experimentally validate the performance of the developed method. Phantoms were prepared by mixing 40 mL of 20% Intralipid (IL) (Sigma-Aldrich; MO, USA) with 750 mL of de-ionized water while absorption was independently introduced by serial additions of 250 mg of dry bovine hemoglobin (Hb) (H3760; Sigma-Aldrich; MO, USA) to the above solution. Liquid phantoms were prepared in a cylindrical container (radius = 6 cm, height = 12 cm) and preliminary experiments showed finite-boundary effects at low scattering solutions (< ∼1cm^-1) and shallow depths (< ∼4 cm). Therefore, the minimum µ_s’ used for the experiments was 5 cm^-1 and the minimum distance between optical channel and container wall was ensured to be greater than 5 cm to best approximate a semi-infinite medium as modeled by the MC lookup tables.

Phantoms were created to have constant scattering coefficient with ∼5% IL by volume [28,29]. Since the optical properties of IL have been well described to be stable and reproducible [22,23], we could also use reference values from previous reports [13,22,23,25] to directly compare the reduced scattering coefficient we estimated to those reported previously. Hb was added five times to create a set of 6 different phantom media having fixed scattering increasing absorption coefficients (corresponding to hemoglobin concentrations ranging between 0-250 µM). This produced phantoms with expected µ_a values varying between 0.003-0.4 cm^-1 for 650 nm, and between 0.022-0.14 cm^-1 for 800 nm.

Each phantom solution was prepared in the container and mixed gently with a magnetic stirrer during measurements. The sensing head of the probe (mounted on a custom probe holder) was lowered until the surface tension of the solution was broken. The probe was clamped in place and the following sequence of measurements were obtained for each phantom. First, the laser wavelength was set to 650 ± 5 nm, the optical fiber for the channel with SDS of 5 mm was manually coupled to the SPAD (via a standard SMA fiber socket) and three repeated scans were acquired. Next, the bandpass was adjusted to illuminate the sample at 700 ± 5 nm, 750 ± 5 nm and 800 ± 5 nm and three repeated TRS measurements were acquired. The same sequence of steps was repeated by manually coupling fibers for the 10 and 15 mm SDS channels into the SPAD, across the same four source illumination wavelengths. Each TRS acquisition was obtained by allowing the TCSPC signal to acquire the signal for 30 s.

Thus, a total of 36 TRS measurements were obtained (3 scans per wavelength, for 4 wavelengths across 3 channels), for each phantom solution. TRS scans collected at 650 nm are shown for all detection channels used (with SDS of 5, 10 and 15 mm) in Fig. 1(b). The difference between t_max of measured DTOF at 5 and 10 mm SDS is also shown and labelled as Δt_{5, 10} in Fig. 1(c). A similar naming convention was adopted when peak differences were calculated using SDS of 5 and 15 mm (Δt_{5, 15}) and SDS of 10 and 15 mm (Δt_{10, 15}), as denoted in Fig. 1(c). Experimental DTOF measurements were obtained with time-resolutions of 12 ps.

3. Results

As described above, the MC lookup table was used to convert the measured Δt (for a given SDS pair) into µ_s’, given an estimate of µ_a. We first tested the developed method to extract µ_s’ for phantom that was mostly scattering with very low absorption – i.e. a solution of 5% IL in water. Figure 3 shows extracted µ_s’ as a function of illumination wavelength for the phantom with 0 µM hemoglobin (5% IL in water) using the lookup method at all three experimental SDS pairs possible (circles: Δt_{5, 10}; diamonds: Δt_{5, 15}; squares: Δt_{10, 15}; stars: literature average from Refs. [13,22,23]). Input values for µ_a were estimated using the tail method from measured data. Derived µ_s’ values in Fig. 3 show expected wavelength dependent decreases in scattering for all the three SDS pairs and tracked expected values from literature. Although Δt_{5, 10} produced µ_s’ closest to reference values, both Δt_{5, 15} and Δt_{10, 15} showed agreement (to better than 25%) with the reported intrinsic scattering parameters of IL [22]. Larger errors are observed at smaller SDS as the uncertainty in derived µ_s’ was high due to the uncertainty in correctly estimating the peak time with the experimental temporal resolution of 12 ps (TCSPC bin-width). Error bars shown represent the standard error from three repeated measurements for each concentration at each wavelength and SDS.

 

Fig. 3. Estimated µ_s’ values obtained from experimental Δt measurements in a solution of 5% IL in water (from each SDS pair) at each wavelength, together with expected (literature) µ_s’ values. Markers show the average of three repeated scans. Data obtained using Δt_{5, 10} most closely matched expected values while those obtained from using Δt_{5, 15} and Δt_{10, 15} still closely followed spectral trends predicted from literature values. The look up estimates were made by obtaining absorption using the tail method (see text). The points are jittered if they share the same independent variable. Error bars represent the standard error for the three repeated measurements at each point (only one side of the error bar is shown for the derived values in the figure for clarity).

Download Full Size | PDF

The values of µ_a used to derive µ_s’ shown in Fig. 3 were estimated using the tail method and could yield inaccurate results. We therefore examined the impact of changing the input µ_a on the extracted value of µ_s’, for any given DTOF. This was done by first, using an initial µ_a0 (obtained for example using the tail method) to derive the scattering coefficient µ_s0’ using the lookup method. This value of µ_a0 was then decreased by 20 - 100% in five equal intervals (where a decrease of 100% was the smallest simulated µ_a = 1.0×10⁻³ cm^-1) to determine an updated absorption coefficient µ_a1 and then using µ_a1 as the input to the lookup method to determine the corresponding updated scattering coefficient µ_s1’.

Percent errors in µ_s1’ (relative to µ_s0’) are shown in Fig. 4, as a function of percent changes in input µ_a0, for the three experimental SDS pairs (bars) used here. Data in Fig. 4 represent DTOFs obtained from four different (representative) media – Medium 1 (Fig. 4(a)): µ_s0’ = 10, µ_a0 = 0.05 cm^-1; Medium 2 (Fig. 4(b)): µ_s0’ = 10, µ_a0 = 0.16 cm^-1, Medium 3 (Fig. 4(c)): µ_s0’ = 13, µ_a0 = 0.05 cm^-1; Medium 4 (Fig. 4(d)): µ_s0’ = 13, µ_a0 = 0.16 cm^-1). DTOFs to estimate µ_s0’ in Fig. 4(a) and 4(c) were obtained from the pure scattering phantom (no hemoglobin) for illumination with 800 and 650 nm, respectively. DTOF used in Fig. 4(b) corresponds to the phantom with 256 µM Hb for 800 nm illumination, while DTOF for Medium 4 (Fig. 4(d)) was from the phantom with 53 µM Hb with illumination at 650 nm (these phantoms were selected because their absorption coefficients for the wavelengths and concentrations used were comparable to each other).

 

Fig. 4. The change in the estimated scattering coefficient given a decrease in the initial absorption coefficient used in the lookup table is shown. Four different media are shown representing optical properties at 800 nm (µ_s’ = 10 cm^-1; (a) and (b)) and 650 nm (µ_s’ = 13 cm^-1; (c) and (d)). µ_a0 values used as input into the lookup table was 0.05 cm^-1 for (a) and (c) while µ_a = 0.16 cm^-1 for (b) and (d). Input absorption coefficient were decreased by 20% of the initial value in five steps and the percent change in retrieved µ_s’ are shown.

Download Full Size | PDF

It is clear to see from Fig. 4 that the derived µ_s’ using Δt_{5, 10} varied less than 8% as the input µ_a values were varied by 100% in media with lower absorption (Fig. 4(a) and 4(c)), as was expected from the MC simulations. Additionally, even for the longest SDS pair (10-15 mm) used, the change in extracted µ_s’ was lower than 20% as the input absorption coefficients changed by 100%. However, when µ_s0’ was obtained with higher absorption coefficients (µ_a0 of 0.16 cm^-1, Fig. 4(b) and 4(d)) the changes in the medium’s absorption translated to larger variations in extracted µ_s’. In general, the larger the error (or uncertainty) in knowledge of a medium’s true µ_a the larger the error in the recovered µ_s’.

Next, we tested the performance of the lookup technique to estimate µ_s’ in phantoms as the absorption coefficient was varied. Acquired DTOF’s using all available SDS and illumination wavelengths (with three SDS pairs and four wavelengths) were analyzed using the MC lookup method to translate the measured values of Δt into µ_s’, in each of the 6 experimental phantoms prepared. Figure 5(a) shows the µ_s’ extracted by the lookup method (for each SDS pair) in phantoms with varying hemoglobin concentrations for illumination with 750 nm and obtaining the input µ_a coefficients from the DTOFs using the tail method. Figure 5(b) shows these data for the same phantoms shown in Fig. 5(a) but used the transmittance method (described in Section 2.4) to determine input µ_a. Both Fig. 5(a) and 5(b) indicate that that the µ_s’ values obtained from Δt_5,10 with the lookup method matched the expected (literature) values (shown as dashed horizontal lines [13,22,23]). Further, all three SDS pairs extracted consistent µ_s’ values in all the six phantom media that had different absorption properties.

 

Fig. 5. Estimated µ_s’ at 750 nm are shown across six concentrations of hemoglobin utilizing two separate methods to determine the µ_a used in the lookup table: (a) the slope of the DTOF tail and (b) measured transmittance of pure hemoglobin. The true (expected) µ_s’ is shown as a dashed horizontal black line. Estimated values using Δt_{5, 10} most closely tracked expected values, but all three SDS pairs once again produced consistent estimates of µ_s’ across the six phantoms. Points are shown as jittered if they share the same independent variable and only one side of error bar are shown here for clarity. Error bars represent the standard error of the estimated scattering coefficient by using each of the three repeated scans for each SDS (in the SDS pair) to estimate Δt.

Download Full Size | PDF

Table 1 shows the mean percent errors between lookup table derived µ_s’ and literature values averaged across all wavelengths, for each SDS pair while using the tail method to determine the µ_a used in the lookup table. Percent errors reported in Table 1 varied lesser than 7% of the values obtained when the transmittance method was used to estimate µ_a. Largest deviations between the two approaches were observed in media with highest absorption where the influence of the IRF impacts the derived values most significantly.

Table 1. Percent errors between estimated and expected (true) values of µ_s’ for each source-detector pair averaged across all experimentally measured wavelengths. µ_s’ was recovered using the decay rate of the DTOF tail to determine the µ_a used in the MC lookup table (see text).

View Table

At wavelengths of 750 and 800 nm, estimating µ_a via the tail method provided improved accuracy in estimating µ_s’. However, at wavelengths of 650 and 700 nm, calculating µ_a using the tail method decreased accuracy in estimating µ_s’ (µ_a obtained via transmittance data yielded lower errors in extracted µ_s’ in these cases). Best accuracy was obtained for media with low µ_a and by using SDS pairs closest to the source (these data cells are highlighted in Table 1).

4. Discussion

We have described a heuristic technique capable of assessing the reduced scattering coefficient of a homogenous medium directly from experimentally measured TRS signals. We have shown that an experimentally measured peak-time difference between DTOFs at two different SDS can be translated to obtain a robust estimate of the medium’s scattering coefficient using preconstructed MC lookup tables. Although our method requires an input (estimated) value for the absorption coefficient, the derived µ_s’ using our approach is only weakly dependent on the input absorption coefficient. Application to experiments with tissue phantoms showed that we could estimate the expected (true) scattering coefficients in media with errors ranging from 5% -25%.

A principal advantage of the presented technique is its inherent speed and simplicity of use (i.e. directly being able to translate measurements into the reduced scattering coefficient of the medium). However, it is worth noting that there are important instrumentation calibration and prerequisite conditions that must be satisfied for it to function accurately. Although our method does not use the IRF for translating measurements into the medium's scattering coefficient, it is not independent of it. It is critical for the IRF (i.e. the intrinsic temporal shape of the incident laser pulse) to remain stable across the duration of experiments and also be identical for each measurement channel used. In our case, each channel used the same length and type of fiber and were all detected by a common detector and as shown in Fig. 1(c) were identical across channels for all wavelengths. Systems that use different fiber lengths and/or employ multiple photodetectors that seek to use this technique might need to consider such issues carefully during development of lookup tables.

Phantom results showed that errors between literature (expected) values and derived values were (on average) lower than 6% when data from SDS of 5 and 10 mm were used and were as high as 25% when data from SDS of 10 and 15 mm were used. Thus, utilizing smaller SDS would yield the best estimates of scattering, consistent with previous reports [5,30–32]. We note that accurately resolving peak time-differences for small SDS requires higher temporal resolution of measured data. Further, non-linear distortions in t_max (from the convolved IRF response) could impact longer SDS measurements more and may explain our findings.

Our findings are consistent with previous reports showing that the peak arrival time (t_max) of the DTOF is strongly sensitive to µ_s’ and only weakly to µ_a [15,17], especially for shorter SDS. It also builds off separate work that seek to mitigate the influence of the IRF in TRS analysis [8,11,33]. Data shown in Fig. 2(a) and 2(b) reflect these findings, where both µ_s’ and µ_a impact the peak time to varying degrees, but as denoted in Fig. 2(c) (shaded bar), a fixed value for Δt (for a given SDS pair) tightly binds the range of allowed µ_s’ values, but not µ_a. Although the technique needs an input estimate for µ_a, this can be derived by applying Eq. (1) directly to the measured DTOF. Accurately using Eq. (1) to estimate µ_a would require a large dynamic range in the detecting system and deconvolving the IRF.

Here, we have shown (Fig. 4) that only a rough estimate of µ_a is needed as it weakly affects the peak position. Utilizing smaller SDS decreases the effect of absorption on the peak time difference leading to more accurate estimates of µ_s’ by reducing the probed volume. Therefore, small SDS would be more sensitive to shallower layers of the medium, while longer SDS could differentially be used to probe deeper tissues. The extension of this method to inhomogeneous media will be addressed in future investigations. Although not experimentally studied here, any fiber geometry with at least two detector channels with unique SDS would be amenable for use with the lookup method presented. Further, using this method to derive the scattering by setting µ_a=0 could be used to obtain a lower bound on µ_s’.

We also note that the SDS must be well defined (within 0.5 mm) to accurately estimate µ_s’. Peak times between DTOF’s at 10 and 11 mm SDS could vary by ∼ 50% translating into recovered error of ∼ 80% in µ_s’. Additionally, although MC simulations and experimental geometries must be matched (i.e. we need experimental measurements from large volumes to represent semi-infinite media) we note a specific advantage of this method in “finite” geometries as only the earliest arriving photons determine Δt – and thus, could be robust boundary effects.

The described method could be used to bootstrap (or confine the range of µ_s’) other analytical inverse solvers to quantify optical properties from TRS measurements, especially in weakly scattering media where DT-based approaches are typically more vulnerable to higher errors [1–3,15,17]. Uncertainty in the index of refraction of the medium also weakly affects the relative peak-time difference. Lastly, since this method works by establishing the temporal peak of measured DTOF reflectance profiles, it should be possible to sparsely sample the DTOF histogram to acquire data faster. These topics will be subject of future investigations.

5. Conclusion

An easy-to-use approach to extract µ_s’ of turbid media that overcomes current limitations in TRS analysis was presented. By using a pair of short SDS (<10 mm), we show that µ_s’ can be estimated within 6% of reference values with reasonable estimations of µ_a. The approach could be further extended to include timing differences of later arriving photons as well as aid traditional inverse solvers in further parameterizing optical properties. The method shows great promise in recovering a medium’s optical transport properties via TRS in regimes where common techniques fail, while facilitating quantitation in real-time without directly measuring the IRF.