1. Introduction

Visible and near infrared diffuse optical imaging (DOI) and spectroscopy (DOS) are used to non-invasively study biological tissues in laboratory and clinical applications [1,2]. Typical tissues investigated include, but are not limited to, brain, skeletal muscles, abdomen, lung, breast and thyroid [3,4]. Measurements are often performed in reflectance geometry, taking advantage of light scattering in highly diffusive media. The reflectance configuration allows using a couple of optic fibers, one for injection and one for light collection, placed on the same surface of the tissue at a distance of few centimeters. The optical properties of the biological tissues, namely the absorption coefficient (µ_a) and the reduced scattering coefficient (µ′_s), can then be estimated with physical models derived from the Radiative Transfer Theory under the Diffusion approximation [5].

In the time resolved approach to DOI and DOS, pulsed laser sources, with a pulse duration of hundreds of picoseconds, are utilized together with fast detectors and timing electronics to record the distribution of the photons time-of-flight (DTOF) when they travel through biological tissues. The use of the time resolved approach, as compared to the widely spread steady state or continuous wave (CW) approach, results in better decoupling when estimating the absorption coefficient (µ_a) and the reduced scattering coefficient (µ′_s), and leads to a straightforward link between photon arrival time and the depth probed by photons [6,7]. The incessant advances in performances of time domain (TD) components and systems, that pave the way to miniaturization and scalability of TD set-ups [8], make this approach the most promising for future applications of DOI and DOS.

A critical step in the TD approach is the need to record the so-called instrument response function (IRF), that is the DTOF when the sample is removed and the injection and collection fibers are directly faced. The IRF is crucial for evaluating of the correct timing of the DTOF, and the IRF characteristics (e.g. width and stability) affect the overall performances of a TD system [9]. In an ideal system, the temporal width of the IRF should be kept as narrow as possible in order to approximate a Dirac delta function. However, in real time resolved systems, where different tradeoffs have to be taken into account, there are many effects that contribute to the IRF temporal broadening, e.g. the finite duration of laser pulses, the temporal dispersion in optical systems, the jitter time of detectors and acquisition electronics. For these reasons, in the interpretation of the TD measurements, the effect of the IRF is usually taken into account by convolution or deconvolution strategies [10,11].

To record an IRF, injection and detection fibers should be placed in a configuration where they are one in front of the other and in contact: we will call this the “ideal geometry” (see Fig. 1(a)). The temporal position of the IRF (usually estimated as the peak or the barycenter of the IRF) is taken as the time origin for the DTOF (i.e. T₀ = 0 ps). This ideal geometry, however, is often replaced by a more complex set-up comprising an attenuation stage and a diffuser (i.e. a thin diffusive medium like a paper foil or a Teflon layer) inserted between the injection and the collection fibers: we will call it the “reference geometry” (see Fig. 1(b)). This configuration is used to protect the detector from excessive incident power and to illuminate all the waveguide modes of the collection fiber (typically a multimode fiber or fiber bundle), mimicking light collection from a diffusive sample and not from a collimated laser beam [12]. The temporal position of the IRF recorded in the reference geometry is delayed with respect to the ideal one because of the additional distance s travelled by photons. However, this delay can be often neglected, due to the very small thickness of attenuation and diffuse layers (usually few tens of micrometer), or easily estimated considering the speed of light in the attenuation material.

 

Fig. 1. IRF acquisition geometries: a) “ideal”: injection and detection fibers are placed one in front of the other; b) “reference”: fibers are placed at a distance s with extra layers in the middle for attenuation and diffusion purpose.

Download Full Size | PDF

The IRF should be acquired for each measurement session in order to be able to monitor possible time drifts that are detrimental for the estimation of the optical properties [13]. However, this procedure can be time consuming and not feasible, especially when performing measurements in clinical environment. Furthermore, there might be no possibility to face injection and collection fibers like in the reference geometry, since they are often blocked (or glued) into a probe designed for measurements in a reflectance scheme [14].

To overcome these problems, we investigated the possibility of recording the IRF in a reflectance scheme without the need to remove the optical fibers from the probe. In section 2 we propose a reflectance configuration for acquisition of the IRF from a reflective surface. Besides, we introduce two different methods to estimate the temporal delay T₁ that occurs when injection and collection fibers are not facing each other to properly set the time scale of the DTOF. Then, in section 3, we validate the proposed methodology through measurements on calibrated phantoms. Finally, in section 4 we critically discuss the findings.

2. Materials and methods

A series of preliminary and validation measurements performed to assess the reliability of the reflectance geometry approach with respect to the reference one are here illustrated. Moreover, two different approaches to determine the delay introduced in the reflectance geometry to be taken into account in TD data analysis are here proposed.

2.1 “Reference geometry” IRF acquisition system

The IRFs in the reference geometry have been acquired with the TD near infrared spectroscopy (NIRS) system described by Re et al. [15]. This device allows the injection of light at two different wavelengths (690 nm and 830 nm) through a multimode glass graded index optical fiber (core/cladding of 100/140 µm and 0.22 numerical aperture (NA)). The injected beam is reflected at 90° by means of a glass prism, obtaining an injection spot of 5 mm diameter. The collection consists of two glass multimode step index optical fiber bundles (bundle diameter: 3 mm, NA = 0.57), with 90° bended terminations. The instrument is equipped with a custom probe able to host one injection fiber and two 90° bended collection bundles at two different source-detector distances (ρ₁ = 15 mm; ρ₂ = 30 mm). The 3D printed probe is presented in Fig. 2(a), while in Fig. 2(b) the picture of the IRF acquisition in reference geometry is shown.

 

Fig. 2. (a) 3D printed custom probe of the TD-NIRS device able to hold one injection and two detection optical fibers. (b) Implemented system for “reference” IRF acquisition and a schematic view of the ray propagation (red arrows) inside it.

Download Full Size | PDF

2.2 Apparatus and methods for IRF acquisition in reflectance geometry

A schematic representation of the novel apparatus for the IRF acquisition in reflectance geometry is reported in Fig. 3(a). As in the reference geometry, stages for attenuation and diffusion are used, trapped between the fibers probe and the frame of the IRF holder. The reflective surface is obtained placing a corrugated aluminum foil at a distance d from injection and collection fibers plane. The aluminum foil was chosen because it is a commonly used reflector, easy to find and to be shaped accordingly to the desired geometry. As stated in a work of 2008 by Janecek and Moses [16], the aluminum foil exhibits specular reflection properties. However, considering the limited range of angles of light emission/acceptance of optical fibers usually employed in a TD NIRS system, an ideal mirror-like surface would limit the number of photons that can reach the detection bundle tips at interfiber distances in the order of a few centimeters. For this reason, we corrugated the aluminum foil to increase the reflection angles and detected photons.

 

Fig. 3. (a) Apparatus for the IRF acquisition in “reflectance geometry”: front view scheme. Injection and collection fibers are on the same plane and rays (red arrows) are collected after being reflected by a corrugated reflective surface at distance d. ρ₁ and ρ₂ are the source-detector separations. (b) “Reflectance” IRF 3D printed holder custom designed for probe in Fig. 2(a). An aluminum foil at the bottom of the box is used as reflective surface.

Download Full Size | PDF

Moreover, in this configuration, a temporal shift T₁ needs to be taken into account to properly set the time scale of the DTOF. The delay T₁ is mainly due to the distance d between the fiber plane and the reflective surface. Nonetheless, the employment of the corrugated aluminum foil to increase the number of reflection angles does not allow to easily calculate the path followed by the detected photons by means of geometrical considerations. Hence, T₁ cannot be directly obtained with a simple distance measurement, but it needs a specific calibration.

Like the “reference” IRF geometry, also the “reflectance” IRF one uses a Teflon layer in front of the detection bundles to excite all the propagation mode of the collection fibers, in a similar fashion to what happens in diffusive media. Therefore, in this work, two different approaches to retrieve empirically the temporal delay T₁ are proposed and schematically depicted in Fig. 4.

 

Fig. 4. Schematic representation of the two approaches proposed to determine the temporal delay of the IRF measured in “reflectance” geometry. The position of the “reference” IRF barycenter has been chosen as the origin of the DTOFs time axis. (a) First method: the “reflectance” delay T₁ is calculated as the difference between the barycenter of the “reflectance” IRF and the barycenter of the “reference” IRF. (b) Second method: µ_a and µ′_s of a calibrated phantom are known and the time shift T₂ is estimated through the DTOF fitting procedure.

Download Full Size | PDF

In the first method, the acquisition of an IRF in reference geometry allows to estimate the delay T₁. Indeed, it can be calculated by estimating the difference between the barycenter of the IRF in reflectance geometry and the barycenter of the IRF in reference geometry (the latter usually assumed as the origin of times for the DTOF).

When it is not possible to acquire an IRF in reference geometry, a second method that exploits a priori information about the optical properties of a previously characterized phantom is here proposed. The IRF in the reflectance geometry is convoluted with a model for TD reflectance [17] to fit the DTOF acquired on a calibrated phantom at a given source-detector distance. The temporal position of the IRF is assumed as the only free fitting parameter, while all other parameters are fixed, in particular µ_a and µ′_s, which are settled to their calibrated values. The delay T₂ is therefore equal to the optimal time shift T_fit determined by the fitting procedure.

Finally, a custom compact apparatus for IRF acquisition in reflectance geometry has been developed. The proposed system has been 3D printed with a black PLA filament (3D Italy) by exploiting a fused filament printer (FDM, Sharebot NG, Sharebot S.r.l., Italy). In this way, we were able to create a custom IRF holder (Fig. 3(b)), tailored as a small box (53×36×33 mm), and suited for the probe previously discussed. For the reasons previously explained, the bottom side of the IRF box has been covered with the corrugated aluminum foil to create the reflective surface. To define the appropriate holder dimensions, and in particular the distance between optical fibers and reflective surface (d), a preliminary characterization to investigate the influence of d on IRF’s barycenter position (t_bar), along with preservation of DTOF shape and full width at half maximum (FWHM), has been carried out. After testing different values ranging from 5 to 35 mm with 5 mm steps, the distance d has been set to 25 mm. Detailed results of this first assessment will be reported in paragraph 3.1.

2.3 Estimation of the time delay

A series of measurements to assess and compare the performances of the methods to estimate the time shift T₁ and T₂ have been carried out. The “reference” and “reflectance” IRFs have been acquired with the holders depicted in Fig. 2(b) and Fig. 3(b), respectively. Furthermore, the delay T₂ has been retrieved by fitting the TD curves acquired on a calibration phantom (µ_a = 0.1 cm⁻¹ and µ′_s = 10 cm⁻¹ at 660 nm). For better performance of the fitting procedure, the parameters optimized were the searched T₂ and µ_a, while µ′_s was kept constant. Moreover, to guarantee the correct timing of the DTOF curves, the position of the IRF barycenter was considered as the time origin of the measurement, like normally done in Time Correlated Single Photon Counting (TCSPC) measurements. In fact, even though the time origin is usually assigned to the maximum of the IRF, the temporal position of the curve peak can suffer from instability due to measurement noise, hence the first moment of the distribution was preferred.

The “reflectance” IRF box has been tested for both interfiber distances ρ₁ = 15 mm and ρ₂ = 30 mm that are present in the probe of the used device (see Fig. 2(a)). For each curve, approximately 10⁶ photons have been acquired in an integration time of 1 s at each wavelength (690 and 830 nm). In order to assess the repeatability of the measurements and to exclude any influence of the Teflon positioning, 10 repetitions for both “reference” and “reflectance” IRF have been performed, changing the Teflon layer each time.

For the validation of the T₁ and T₂ estimation approaches, the same protocol designed to test the reflectance geometry IRF acquisition, in terms of number of counts per seconds acquired, integration time, wavelengths and interfiber distances, has been followed to measure a set of homogeneous solid phantoms with known nominal optical properties [18]. Eight phantoms with fixed µ′_s equal to 10 cm⁻¹ and µ_a linearly varying from 0.01 cm⁻¹ to 0.49 cm⁻¹ in steps of 0.07 cm⁻¹ at 660 nm have been tested to assess the influence of T₁ and T₂ over the quantification of the absorption coefficient. On the other hand, in order to investigate the effects of T₁ and T₂ on the estimation of the reduced scattering coefficient, a series of four phantoms with constant nominal µ_a = 0.07 cm⁻¹ and linearly changing µ′_s from 5 cm⁻¹ to 20 cm⁻¹ in steps of 5 cm⁻¹at 660 nm has been measured. The fitting procedure of the measured DTOF acquired on these solid phantoms has been applied using firstly the “reference” IRF, then the “reflectance” IRF with both the T₁ and T₂.

3. Results

In the following section, results obtained for both the preliminary characterization and validation of the reflectance geometry for IRF acquisition in TD measurements are described.

3.1 Preliminary characterization

In order to develop a probe for the acquisition of the IRF in reflectance geometry suitable for a TD NIRS device, a set of preliminary measurement has been carried out. A first test was aimed to assess if the novel acquisition geometry introduces any distortion in the IRF shape and if the FWHM is preserved. In Fig. 5(a), “reflectance” IRF curves acquired at ρ₂ = 30 mm at both wavelengths have been overlapped to the corresponding “reference” one. No distortion of the shape can be observed, although a slight average increase of FWHM of 23 ± 4 (31 ± 7) ps at 690 (830) nm was shown. The same overall good overlapping was guaranteed for the curves acquired at source-detector distance ρ₁ = 15 mm (graph not reported), but a small decrease of FWHM was found, with an average difference of 19 ± 5 (13 ± 4) ps at 690 (830) nm. Moreover, to verify that the chosen distance between fibers and reflective surfaces does not introduce any non-linearity in the photons arrival time, namely in the barycenter position (t_bar) of the IRF, a series of measurements with d ranging from 5 to 35 mm has been performed. In Fig. 5(b), the relationship between t_bar and d is depicted for the two wavelengths (690 and 830 nm) and source-detector distance ρ₂ = 30 mm. As expected, when increasing the distance between the fibers plane and the reflective surface, the arrival time of photons increases. It is important to highlight that the relationship between t_bar and d is linear as proved by the linear interpolation in Fig. 5(b) (R² > 0.99 for both series). Hence, no distortion due to the extra pathway followed by photons is introduced. Comparable results have been obtained for ρ₁ = 15 mm.

 

Fig. 5. (a) Comparison between “reference” IRFs (red → 690 nm, purple → 830 nm) and “reflectance” IRFs (blue → 690 nm, green → 830 nm). (b) IRF barycenters in reflectance geometry as a function of the distance d between fibers and reflective surface at 690 nm and 830 nm and ρ = 30 mm. The two wavelengths are shifted along the time scale in order to avoid overlap and better appreciate them.

Download Full Size | PDF

In conclusion, no particular limitation about the distance between fibers and reflective surface needs to be considered in the design of a “reflectance” IRF holder for a TD device with comparable optical fibers and NA to the one tested here. Considering the overall dimensions, a 3D printed box with a distance d = 25 mm was chosen (see Fig. 3(b)).

3.2 “Reflectance” time shift estimation and validation of the methods

In this paragraph, the results obtained for the two approaches are reported, followed by the results of the validation measurements comparing both “reference” and “reflectance” IRFs performances in the fitting procedure for the estimation of the optical properties. In Table 1, time shift T₁ and T₂ values of the “reflectance” IRF are reported.

Table 1. Average values and standard deviations of T₁ and T₂ calculated with the two methods, at two different source-detector distances ρ and for two wavelengths. Statistics of the values are based on 10 repeated measurements. For T₂ estimation, the µ′_s values set to initialize the fitting procedure was 9.32 cm⁻¹ at 690 nm and 7.35 cm⁻¹ at 830 nm.

View Table

As expected, for both methods, higher values of the time shift have been found for the longer interfiber distance ρ₂. This is consistent with the fact that increasing ρ, the path covered by reflected photons also increases. Additionally, a small dependence of T₁ and T₂ on the wavelength can be observed for the short interfiber distance ρ₁. Moreover, differences between results from the two methods can be observed. However, considering the standard deviation reported in Table 1 and the 3.05 ps/channel resolution of the TD instrument used, we would not expect these discrepancies to greatly affect the retrieved absolute values of µ_a and µ′_s. To check this, a further assessment of the influence of T₁ and T₂ upon the fitting procedure is necessary.

Consequently, the calculated time shifts have been used to calculate the optical properties of the solid phantoms tested as described in paragraph 2.3. In Fig. 6, a comparison of µ_a and µ′_s values estimated using the “reflectance” IRF against the ones obtained with the “reference” IRF is reported. The results shown refer to the longer interfiber distance ρ₂ at both wavelengths and they represent the average values obtained over 10 repeated measurements. The small error bars depicted represent their standard deviations. No substantial difference seems to be introduced in the estimation of µ_a and µ′_s by the usage of the two different time delays (T₁ results in blue and T₂ red diamonds). Besides, the estimated values do not deviate significantly from the superimposed axis bisector obtained using the values estimated with the “reference” IRF, indicating that, for increasing µ_a and µ′_s, the desired linearity in estimation is preserved.

 

Fig. 6. Phantoms optical properties (µ_a in the first row and µ′_s in the second row) at two wavelengths (columns) and ρ₂ = 30 mm. Average values over 10 repeated measurements and their respective standard deviations (error bars) are reported. Optical properties estimated using a “reflectance” IRF corrected with two different calculated time shifts (T₁ in blue and T₂ in red) against the ones obtained with the “reference” IRF (solid black line) are plotted.

Download Full Size | PDF

Finally, we calculated the relative error of the optical properties, ɛ_µa and ɛ_µs, of “reflectance” approaches with respect to the “reference” one as:

In Fig. 7, these relative errors are reported as a function of nominal absorption coefficient values. The results presented here refer only to phantoms with the same µ′_s equal to 10 cm⁻¹. As can be observed in the first column of Fig. 7, an error lower than 10% is always obtained for both µ_a or µ′_s when determined by fitting the phantoms DTOF using T₁. It is also noteworthy that the accuracy decreases for higher absorption values and for the shorter interfiber distance. As an example, with a 3 cm source-detector separation, the ɛ_µa (ɛ_µs) is reduced to 3% (4%).

 

Fig. 7. Relative errors for µ_a (first row) and µ′_s (second row) obtained with a “reflectance” IRF with respect to the “reference” one as a function of the nominal absorption coefficients (constant µ′_s = 10 cm⁻¹) . Both wavelengths and source-detector separations are considered. On the left hand side, the graphs show the results obtained with the time shift T₁, while the ones on the right hand side show results relative to T₂.

Download Full Size | PDF

Comparable results have been achieved for the accuracy of coefficients estimated by means of T₂, even if the reduced scattering coefficient appears to be less affected by the method (ɛ_µs < 4% at 690 nm and < 2% at 830 nm).

Concerning the estimation of optical properties of phantoms with constant µ_a equal to 0.07 cm⁻¹ and different µ′_s (5, 10, 15, 20 cm⁻¹), relative errors below 2% for µ_a with T₁ and T₂ has been found for both interfiber distances and wavelengths (data not shown). The second method appears to be less accurate in the estimation of µ′_s, especially for lower scattering values, i.e. 5 cm⁻¹, at short interfiber distance ρ₁ = 15 mm. In these instances, the obtained ɛ_µs are < 9.5%, while in the other cases ɛ_µs is lower the 2%.

4. Discussion and conclusions

In TD DOI and DOS, the acquisition of the IRF is crucial for the correct timing of the DTOF and for the overall performance of the system. The introduction of a reflectance geometry for the acquisition of IRF would be greatly beneficial in all those cases where a probe realized for reflectance measurements can not be modified without increasing the complexity of the measurement protocol. To the best of our knowledge, only for one hybrid Diffuse Correlation Spectroscopy/Time Resolved NIRS device described in the work published also by some coauthors of this work (Giovannella et al. [14]) a similar approach for IRF measurement in reflectance geometry has been used. However, we are not aware of the presence in literature of a systematic validation of this method, demonstrating that it can be applied without introducing errors into the estimation of the medium optical properties. The aim of this paper is to demonstrate the reliability of this approach, while making the TD NIRS users aware of the criticism that can be met during its implementation, e.g. the type of reflective surface needed according to the optical fibers used, the influence of the source-detector separations, the estimation of the time delays introduced in the measurement.

With this purpose, considering the nowadays spread of 3D printing materials for different applications and their versatility, which can be useful also for custom-made fiber holder, a 3D printed solution is here presented and characterized. Moreover, particular attention has been devoted in this work to provide the users with two different validated methods to retrieve the right effective time shift introduced by adopting this approach for IRF measurement.

The first important finding that has to be pointed out from the results obtained is that the “reflectance” IRF approach is reliable, provided that a preliminary characterization is performed, and it is feasible to be implemented in a compact system, easy to handle also in a clinical environment from non-expert operators.

The first measurements performed by changing the distance between optical fibers plane and the reflective surface allowed us to assess that no distortion is introduced due to the geometry of the probe at the source-detector separation distances and optical fibers NA tested here. Therefore, a 3D printed holder could be designed for each TD NIRS device and employed in the subsequent measurements. Further investigation demonstrated a good accuracy of this approach in the quantification of the optical properties. The small errors affecting µ_a and µ′_s (ɛ < 10%) found when the “reflectance” IRF is considered instead of the “reference” one, demonstrate that the variation, about 20 ps, registered in the FWHM of the “reflectance” IRF does not affect these estimations significantly.

Moreover, the issue of the estimation of the time shift introduced in the IRF acquisition in reflectance geometry has been addressed with two different approaches. The first approach is based on the availability of a “reference” IRF, while the second one relies on the knowledge of the optical properties of a calibration phantom. The two methods lead to time shift values, specific for wavelength and ρ, that differ less than 10 ps. Provided the temporal resolution of the TD instrument, this error can introduce a negligible misalignment. As reported for the results of the validation measurements over the estimation of µ_a and µ′_s of a set of characterized homogenous solid phantoms, there is not a significant difference between the results obtained with the two approaches. Both of them suffer from higher quantification errors in media with higher absorption and lower scattering. Overall, we found that the estimation accuracy increases for higher source-detector separation, with relative errors compatible with the results obtained in previous works. However, in contrast to findings in literature, the µ′_s estimation seems to be less affected by time uncertainties [19].

Limitations of this study are the use of a set-up with a relatively broad IRF (around 600 ps) and a limited number of wavelengths (690 nm and 830 nm). Indeed, the work was specific for TD NIRS applications targeting brain and/or muscle. Work is in progress to extend the characterization to the case of narrower IRF and broadband or multi-wavelength TD DOI and DOS system.

In conclusion, both the approaches to retrieve the time shift to re-align “reflectance” IRF provide good results in the quantification of optical parameters. If during the preliminary characterization measurements of the selected device the acquisition of a “reference” IRF is feasible at least one time in conjunction with a “reflectance” IRF, the first approach and the calculation of the here so called T₁ as the difference between the barycenters position is suggested. It is sufficient to perform this measurement only once because, even if a drift of the system would occur in subsequent sessions, there is no reason for a change of T₁. For particular circumstances, especially for commercial devices, where the probe is a closed box and it is not possible to unmount the optodes to measure a “reference” IRF, the second method and the estimation of T₂, is equally applicable. In the latter case, a calibration phantom with known optical properties needs to be provided in order to avoid quantification errors. It should be noted that in order to achieve the best performances and minimize time uncertainty errors, the characterization of the shift value has to be specific for every configuration of optical probe (i.e. source-detector distance) and for fiber characteristics (i.e. fiber NA). Once the calibration procedure is carried out, the values obtained revealed to be stable and applicable for measurements acquired in time for that specific device.