1. Introduction

Quantifying the optical properties through biological tissues is especially important for clinical diagnosis [1], owing to the pathological and physiological indicative values that the optical properties can reflect and represent for many essential biomarkers such as oxy-/deoxy-hemoglobin and cytochrome-c-oxidase in tissues (e.g. brain, kidney, and breast [2–4]). However, the measurement of the optical properties is also difficult due to the structurally complex random and optically diffusive nature that tissues have in the near-infrared therapeutic window (600 nm < λ < 1000 nm) [5]. Diffuse optics, namely the implementation of near-infrared light propagation through highly diffuse media to probe the optical properties of the media in depth [6], has been studied through the years for its special advantages compared with other techniques regarding the chromophore-selectivity and the non-invasiveness when retrieving the optical properties in tissues.

The optical properties in diffuse optics, which primarily affect the light propagation and interaction within tissues, are optical absorption and scattering. Absorption is of particular biomedical interest since it carries information about tissue chromophores. However, the transport of near-infrared light in tissues is strongly dominated by scattering, and thus the reconstruction of absorption is heavily hampered by the diffuse propagation of photons. The presence of unknown scatter would cause undesirable artefacts for the recovery of absorbers [7]. Moreover, similar scattered optical fields can be generated from different combinations of optical absorption and scattering properties. In other words, retrieving absorption in diffuse media is an ill-conditioned inverse problem [8]. Therefore, it is crucial but difficult in diffuse optics to simultaneously retrieve absorption and scattering with sufficient accuracy.

Many studies have been devoted to this field following various approaches. In typical diffuse optics measurements, light illuminating the surface of diffuse media and propagating through the media would be collected by a detector or detector array. The variants of such basic scenario are abundant. The light source could be pulsed [9] or continuous [10], spatially structured [11] or multispectral [12]. Also, the detected light signal could be in time domain, in frequency domain [13] or continuous wave. However, regarding the data analysis process, most of these techniques quantify absorption and scattering in similar indirect procedures. By fitting a certain appropriate forward model (e.g. Monte-Carlo method [14], diffusion theory with approximations [15,16]) to the measured data, the inverse problem is considered as solved. These techniques would provide different data types which have intrinsically different level of information content. The choice of techniques and of the corresponding data types profoundly influences the accuracy of the retrieved optical properties. Nevertheless, for measurands in all these data types, such as space resolved (SR) amplitude in space domain (SD), time resolved photon distribution in time domain (TD), or demodulation and phase shift in frequency domain (FD), the effects originating from absorption and scattering of the diffuse media are always heavily entangled and cause crosstalk and ambiguity [17,18] in general. Thus, the simultaneous retrieval of absorption and scattering coefficients with high accuracy remains to be a tricky challenge today, especially for the circumstances of high noise and perturbations when measuring in vivo tissues.

On the other hand, the structural complexity of many biological tissues sets the limitation of applying homogenous models despite their success in infant and small animal studies [19,20]. The assumption of tissue homogeneity would impede robust and reproducible reconstruction of optical properties in the human tissues, and the non-negligible depth dependence of the tissue composition has also been reported [21]. Nevertheless, using real detailed anatomical structures for forward models may also not be necessary since it was shown not to improve retrieval performance compared to slab-based layered models during the fitting of measurements [22]. Therefore, layered models may be a reasonable choice. This is also motivated by the fact that structures of many tissues are layered, for instance muscle tissue underneath superficial fat or the layered structure of the human head. Many studies have already developed the techniques such as TD single distance [23], TD multi-distance [22], FD multi-distance [24], and improved diffusion model for TD and FD dual-distance [25], and have shown good retrieval performance for the absolute quantification of optical properties based on layered models. Another study [26] considered both TD and SD information to fast calculate optical absorption in layered models with optical scattering as prior knowledge. In a previous work [27], we have introduced the first attempt of applying the Space-enhanced Time Domain (SeTD) method in homogenous cases and demonstrated its good performance. The aim of the present paper is to evolve the SeTD method and investigate its feasibility for more practical and clinically relevant layered structures.

In this study, we propose a novel methodology that merges the information acquired from different domains to obtain an optimal estimate of the optical properties in both layers, by applying the SeTD method in two-layered structures. Investigations on two-layered liquid phantoms and analyses by TD, SD, and SeTD methods demonstrate the retrieval performance of these methods regarding accuracy, stability and linearity of optical properties. According to a finding in TD, we introduce the term of “deep scattering neutrality”, which is briefly discussed regarding the reason why it is nearly impossible to recover scattering in the deep layer by TD alone and how the SeTD method can accurately determine it.

2. Concept of space-enhanced time-domain method

With the combination of the information acquired in time and space domains, optical absorption and reduced scattering coefficients (μ_a and $\mu_{\mathrm{s}}^{\prime}$) can be quantified more accurately by using the intrinsic independence of the two properties. A comprehensive explanation of the SeTD approach and its advantages have already been given in a previous article [27] that considered homogenous media. In layered structures the general principle remains the same. Namely, the effects from optical absorption and scattering of diffuse media on the diffuse reflectance are entangled and correlated in time domain (TD) and space domain (SD) in different ways. By the effective mutual combination of TD and SD information, the simultaneous retrieval of absorption and scattering properties can be achieved with improved accuracy, compared to single-domain methods.

In TD, the effects of absorption and scattering on the temporal profiles of distributions of times of flight (DTOF) of photons are opposite to each other, i.e. increasing scattering broadens a DTOF, and increasing absorption narrows it. In order to maintain a certain DTOF profile, the changes of μ_a and $\mu_{\mathrm{s}}^{\prime}$ should occur in the same direction, i.e. there is a positive correlation of them in {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space.

In SD, the shapes of spatially-resolved amplitude curves (SRACs) are mainly determined by the effective attenuation coefficient μ_eff, where μ_eff = (3μ_a $\mu_{\mathrm{s}}^{\prime}$)^1/2. To maintain a certain SRAC shape, for $\mu_{\mathrm{s}}^{\prime}$ and μ_a any increase of one coefficient should be compensated by a decrease of the other, by keeping μ_eff constant [28]. Thus, $\mu_{\mathrm{s}}^{\prime}$ and μ_a follow an approximately reciprocal function under constant μ_eff, i.e. a negative correlation in {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space.

For layered structures, if one focuses on the optical properties in the deep layer and their correlations in TD and SD, the aforementioned general trends and correlations would not change, but only vary to slight extent. For illustration, Monte-Carlo simulations in turbid media of two-layered slab geometry were performed to demonstrate the correlations of scattering and absorption in the second layer and validate the SeTD method. Following the well-known concepts of MC simulations in diffuse optics [29–31], a GPU-based implementation [32] was applied with time-resolved recording of the simulated photons in reflectance. To generate DTOFs, 42.9 billion photons were launched into the two-layered slab where the 1st layer was 10 mm thick and the 2nd layer was 47.8 mm thick. For the 1st layer the optical properties are constant at (μ_a1, $\mu_{\mathrm{s} 1}^{\prime}$) = (0.10, 10.65) cm⁻¹, while for the 2nd layer the optical properties cover the ranges of μ_a2 ∈ [0.05, 0.25] cm⁻¹, $\mu_{\mathrm{s} 2}^{\prime}$ ∈ [6, 14] cm⁻¹.

The MC simulations of a combination (μ_a2, $\mu_{\mathrm{s} 2}^{\prime}$) = (0.15, 10) cm⁻¹ were selected as virtual measurements. Poisson noise was added to the DTOFs at selected detector positions ρ to imitate TD and SD measurements. The total photon count and time bin width were set to 10 million and 2.056 ps, respectively. To derive the optical properties by fitting measurements with simulations, the error norm was applied in χ² form as the universal objective function:

 

Fig. 1. Determining (μ_a, $\mu_{\mathrm{s}}^{\prime}$) of the 2nd layer based on MC simulations. The χ² distributions are plotted for (a) TD, (b) SD, and (c) SeTD methods in {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space. White crosses mark the true values. Results from the methods are marked by green, yellow and magenta crosses, i.e., where the χ² reach the minimum. Red lines are contours where χ² are 3 times of the minimum. The step sizes of the maps are 0.0005 cm⁻¹ for μ_a axis and 0.02 cm⁻¹ for $\mu_{\mathrm{s}}^{\prime}$ axis.

Download Full Size | PDF

From Figs. 1(a) and 1(b) it can be also seen that the discrepancies between the results from TD and SD methods and the true values are remarkable, indicating that neither the TD nor the SD method alone are sufficient to accurately retrieve the deep layer’s optical properties. In contrast, their effective combination could be promising to improve the retrieval accuracy. The SeTD method is to decouple absorption and scattering and reduce the ambiguous area among them in a new artificial spatio-temporal domain (ST). In the ST domain, one optical parameter’s artefacts which are introduced by errors in the estimation of the other parameter, can be minimized. χ_ST², the ST error norm distribution in the {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space, is introduced as:

During the process of iteration, the correlation coefficient of the two parameters (μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$) taken at the group of points located on the contour line of χ_ST² is considered as the objective target criterion. And the iteration optimizes and reduces it as close as possible to zero. If the correlation is almost zero, the retrieved μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ are then nearly independent. By adjusting λ along with iteration, the method is balancing the weight of the information from TD and SD to find the optimal λ for the artificial spatio-temporal domain, whereupon the two unknown parameters can be decoupled. As presented in Fig. 1(c), the effective integration of TD and SD information leads to a more well-defined minimum, i.e., the optimal convergence in {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space. For the case shown in the figure, the retrieved values from SeTD method are (μ_a2, $\mu_{\mathrm{s} 2}^{\prime}$) = (0.152, 9.64) cm⁻¹, which are rather close to the true values.

The orthogonality of both coordinates on the optimized contours implies that the crosstalk among μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ is reduced. The smoothness within the TD contour is overcome by the steepness of the SD contour, indicating that the ambiguity of solutions is reduced. Under the same noise level, the shrinking contour reveals a more confined and stable solution of the inverse problem.

3. Methods and materials

3.1 Experimental setup

The experiments in this study were performed with a laboratory setup, as schematically depicted in Fig. 2. The light source was a supercontinuum laser (SuperK Fianium FIU-15 PP, NKT Photonics, Germany), providing picosecond light pulses with a repetition rate of 39 MHz. A wavelength band centered at 830 nm was selected from the broadband beam by a variable bandpass filter (Varia, NKT Photonics, Germany) with a bandwidth of 10 nm. The near-infrared beam was then collimated into and guided by a 1.6 m multimode graded index fiber (GI-Fiber 1; core diameter 400 μm, NA 0.27, Leoni Fiber Optics, Germany) to the phantom. The photons were scattered and partially absorbed in the phantom. To collect diffusely re-emitted light, another 1.6 m multimode graded index fiber (GI-Fiber 2; core diameter 600 μm, NA 0.22, Leoni Fiber Optics, Germany) was used. The measurements were carried out in reflection geometry, i.e., both fibers’ ends equipped with SMA ferrules were plugged into the holder and vertical to the surface of the phantom cell. The source-detector distance ρ was precisely defined by the positions of cylindrical holes in the NC-machined holders. Both fibers’ ends were carefully positioned in contact with the phantom’s surfaces.

 

Fig. 2. Schematic diagram of the experimental setup. The solid red lines are optical paths and the dotted blue lines are electric signal paths.

Download Full Size | PDF

The detection components were installed in a black metal box to avoid ambient light. The light from GI-Fiber 2 was first collimated; then the free space beam was reshaped by a focusing lens and hit the photocathode of the detector, a hybrid photomultiplier module (HPM-100-50, Becker&Hickl, Germany) operated by a detector control card (DCC-100, Becker&Hickl, Germany). Two shutters (Melles Griot, Germany) and continuously variable neutral-density filters (NDC-100C-4, optical density 0 to 4.0, Thorlabs) were placed in the free space beams at laser and detector, to adapt the light intensity to the appropriate level for detection. Filters were motorized and controlled by a dedicated LabVIEW program.

The photon pulses collected by the detector were recorded by a time-correlated single photon counting (TCSPC) module (SPC-150, Becker&Hickl). The Sync signal from the laser system was transferred to the TCSPC module to provide the timing for the CFD (constant-fraction discriminator) signal from the detector. In this work, DTOFs were recorded with 1 s collection time and 20 repetitions. The dead time compensation mode of TCSPC was activated to prevent potential counting loss due to dead time effects, in case of higher photon exposure which could be inevitable for the measurements at short source-detector distances. Pulse shape distortions due to dead time effects were avoided by limiting the maximum count rate to <10⁶/s. TCSPC control, and real-time display of DTOFs and raw analysis of several characteristics were all performed by another dedicated LabVIEW program.

In addition to DTOF measurements on phantoms, it is crucial to have an independent and accurate instrument response function (IRF) measurement for TD analysis. A cage-system based “response box” was prepared, in which both fibers’ ends are mounted to face each other with a distance of 54 mm. In front of each fiber’s end a sheet of paper is inserted to engage the full numerical aperture of both fibers and to mimic a similar temporal dispersion in the fibers as in the DTOF measurements on phantoms. The optical pathlengths between both fiber’s ends is precisely evaluated to determine t₀, i.e., the time shift between DTOF and IRF measurements. For this laboratory setup, the overall time resolution, represented by the full width at half maximum (FWHM) of the IRF, is about 105 ps at 830 nm.

3.2 Preparation and characterization of the layered phantom

Liquid mixtures of intralipid, ink and water are used as phantom material. Intralipid-20% (Fresenius, Austria) served as diffusive medium for its stability over long time, optical similarity with tissues after proper dilution, and minor variation among batches [33]. Indian ink (Higgins ink #44201, Chartpak, USA) served as absorber, since it is chemically well-dispersible and stable, and practically does not alter the scattering properties of intralipid at the concentrations used [34]. Many researches have investigated these two components’ applicability for diffuse optics experiments [35,36]. It has been shown that a phantom with well-defined optical absorption and scattering properties can be gained through carefully weighing the two components and calculating their concentrations in the mixtures.

By a well-accepted methodology [37], the components were first characterized in terms of absorption coefficient for dilute ink’s concentration, and of reduced scattering coefficient for intralipid’s concentration at the relevant wavelength 830 nm, with independent time-domain measurements. These measurements were performed in transmission and reflection geometries on a homogeneous phantom, i.e. a 3 cm thick black PVC cuboidal cell with small (diameter: 3 mm) transparent circular plexiglass windows [36], containing known concentrations of intralipid and ink. By this characterization, linear relations between concentrations and optical properties are established. Hence, the nominal values of μ_a, and $\mu_{\mathrm{s}}^{\prime}$ of the layered phantom can be estimated from the concentrations during the preparation procedure. In this work, for the layered phantom, the nominal values of first layer’s optical properties (μ_a1, $\mu_{\mathrm{s} 1}^{\prime}$) are fixed at about (0.1, 9.3) cm⁻¹, whereas the second layer’s μ_a2 is varied and investigated in 10 cases in the range [0.03 0.2] cm⁻¹ while $\mu_{\mathrm{s} 2}^{\prime}$ is approximately constant at 9.2∼9.3 cm⁻¹. The slight change on $\mu_{\mathrm{s} 2}^{\prime}$ value is due to the intralipid’s concentration change induced by adding small amount of dilute ink during the measurements. The relative uncertainty of the nominal values is evaluated as approximately 1%. The values for optical properties were chosen to be in a range relevant to brain imaging [38]. Figure 3(b) shows the nominal values of absorption in both layers for the 10 investigated cases. The detailed values are also listed in Table 1 in Section 4.

 

Fig. 3. Layered phantom. (a) PVC cell for phantom; (b) μ_a nominal values of both layers in ten cases; (c) structure of phantom cell and source/detector positions in the measurements.

Download Full Size | PDF

As demonstrated in Fig. 3(c), for two-layered measurements, a 50 mm-thick black PVC cell was prepared as the container of the liquid mixture, where the first layer was separated from the second layer and the PVC front wall by two 50 μm Mylar membranes (DuPont, marked by 2 blue lines). The inner width and height of the entire cell are 120 mm and 135 mm. It has been found that the effect of thin Mylar layers in diffuse reflectance measurements can be considered negligible [39]. On the other hand, due to the flexibility of Mylar membranes, the actual thicknesses of the two layers are influenced by the liquid levels of both layers and must be precisely measured for its importance in the correct forward model. For this purpose, an independent experiment was carried out to assess the time difference of perpendicular light echoes from a laser beam hitting both Mylar membranes when the layers were filled with clear water to certain levels. From these measurements, the actual thicknesses in the two-layer configuration were found to be 11.05 mm for the first and 38.95 mm for the second layer.

A slice of black velvet vinyl is glued on the entire inner surface of the PVC front wall [top in Fig. 3(c), marked by the green line] to prevent stray light travelling along the Mylar surface between injection and collection fibers. The front wall has holes with the same diameter as the fiber ferrules to hold the fibers. Six source-detector distances (10, 25, 30, 35, 40, 60 mm) are included in the diffuse reflectance measurements, for which the assignments are presented in Fig. 3(c). The 10 mm measurements are used for determining the first layer’s properties, while the other measurements are used for determining the second layer’s properties.

3.3 Monte Carlo simulations and forward model

To analyze the measurements on the two-layered phantoms, the GPU-based MC program introduced in Section 2 was applied for forward simulations of the DTOFs remitted at the ρ values used in the experiments. To account for side effects caused by the Mylar foils, a four-layer slab geometry was used in the simulations. Hereby, the first and third layer are the Mylar membranes while the second and fourth layer are the compartments containing dilute intralipid and ink. The thicknesses of the layers are S₁ = 0.05 mm, S₂ = 11.05 mm, S₃ = 0.05 mm, and S₄ = 38.95 mm. The refractive indices of Mylar and intralipid are set as n_m = 1.65 and n_i = 1.33 according to literature [35,39]. The ambient refractive indices are n_a = 1 for air and n_b = 1.54 for the bottom PVC wall. The scattering coefficient of Mylar is assumed as 6 cm⁻¹. The anisotropy factors of Mylar and intralipid are set as g_m = 0.5 and g_i = 0.7. It is worth to note that the g value does not significantly affect the simulated photons’ trajectories for the source-detector distances of interest in this work, whereas the refractive index needs to be assigned correctly since the mismatch at interfaces between different media substantially influences the re-emitted photon’s distributions in both time and space domain. Fresnel reflections are taken into account at all the n-mismatching interfaces. The lateral edges of the phantom cells can be neglected since the measurements were always carried out in the central part of the cell’s frontal surface. The high absorption of PVC and the black velvet vinyl also prevented any uncontrolled light guiding propagations along the surface that can hardly be modelled.

To calculate the DTOFs for a single set of the layered optical properties, 300 billion photons are launched into the medium. Simulated photons re-emerging from the slab’s surface are recorded on concentric rings surrounding the initial light injection point with a step size of 0.1 mm and a time resolution of 4.07 ps (the time channel width in the TCSPC measurements). Then the simulated DTOFs are convolved with the typical Gaussian beam profile according to the diameter of GI-Fiber 1. The cross section and diameter of GI-Fiber 2 are considered as well by overlaying a circular detection area on the simulation output.

A MC simulation database was created in the form of a look-up table (LUT) of regularly chosen (μ_a2, $\mu_{\mathrm{s} 2}^{\prime}$) combinations for the second layer. Hereby, the parameters (μ_a1, $\mu_{\mathrm{s} 1}^{\prime}$) of the first layer were fixed to the values obtained from the measurement at ρ₁ = 10 mm as described below. The range of μ_a2 was [0.02 0.22] cm⁻¹ with a step size of 0.005 cm⁻¹, and the range of $\mu_{\mathrm{s} 2}^{\prime}$ was [4.5 14.5] cm⁻¹ with a step size of 0.25 cm⁻¹. Other possible combinations can be obtained by 2D linear interpolation among the simulated DTOFs of adjacent (μ_a, $\mu_{\mathrm{s}}^{\prime}$) pairs. To separately analyze the measurements for ρ₁ = 10 mm and retrieve the optical properties of the first layer, another MC simulation database was generated by using the geometry of a 50 mm-thick homogeneous slab. These simulations were performed by white MC calculations as described in [40], i.e. with setting μ_a of the slab to 0. The desired absorption values would be calculated outside of the MC simulations then. Convolution for finite laser beam and detector size were also done as described above for the layered situation.

3.4 Measurement method and inverse procedure

In this work a sequential protocol is designed to retrieve the optical properties of the first and the second layers one after the other. The involved steps are as following:

The inverse models for TD and SD were applied by using the χ² objective function in Eq. (1) whereas for the SeTD method the model was based on Eq. (2), to retrieve the optical properties in the second layer. The traversals of all combinations of μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ were carried out to evaluate the discrepancy between measurements and forward models and to draw χ² error norm distributions in {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space. The optimal solutions of all the three models are determined at the place where the objective function becomes minimal.

Measurements at short (ρ₁) and long (ρ₂) distances can be used to distinguish and estimate the optical properties of superficial and deep layers, respectively. This prerequisite was proven by considering the density of the photon trajectories for the different ρ values. It is well-known that this density follows the so-called “banana shape” [41]. The mean trajectories of detected photons can be visually presented as the spatial distribution of the mean partial pathlength (MPP) values. MPP is the mean length that a photon travels inside a defined volume on its way from the source to the detector. Another type of MC simulation [42], which was modeled in accordance with the phantom structure in Section 3.2, was applied to illustrate the trajectories of detected photons at various ρ. This MC simulation in a voxel-based model allows to obtain MPP values for any complex 3D structures with varying optical properties. 5 billion photons were launched for the homogeneous optical properties: μ_a = 0.01 cm⁻¹, $\mu_{\mathrm{s}}^{\prime}$ = 1 cm⁻¹, and n = 1.33. The voxel size was 0.25 mm in each direction. The results for ρ = 10 and 60 mm are presented in Fig. 4. The 2D spatial distribution of MPP values was obtained by summing all voxels perpendicular to the plane locating the source and detector. The distributions demonstrate the expected “banana shape” for short and long distances, i.e., the longer the distance is, the more likely that the detected photons have traveled through the deep region of diffuse media. The 1D depth profiles, obtained by further summing all voxels along the other lateral direction, quantitatively show the depth profiles of MPP values for each 1 mm steps. Almost all photons collected at ρ₁ = 10 mm have travelled in the superficial layer, whereas for ρ₂ = 60 mm a considerable fraction of the collected photons has experienced the region deeper than 10 mm. Therefore, step 1 of the protocol can reliably use the homogenous MC simulations to retrieve the optical properties of the superficial layer, and measurements at step 2 can be used to retrieve the deep layer’s optical properties by using the two-layered MC simulations.

 

Fig. 4. Mean partial pathlength (MPP) for ρ₁ = 10 mm (top row) and ρ₂ = 60 mm (bottom row) of the simulated photons for the layered structure with homogeneous optical properties. (a, c) 2D: MPP values summed along the direction perpendicular to the plane locating the source and detector. (b, d) 1D: Depth profiles of MPP values summed for each 1 mm steps. The purple line highlights the boundary between the first and second layers.

Download Full Size | PDF

4. Results

Table 1. Comparison of nominal values and retrieved values from three approaches (TD: Time Domain, SD: Space Domain, SeTD: Space-enhanced Time Domain) of the second layer’s optical properties in ten absorption-incremental cases. The best matches to the nominal values are marked in bold. (unit: cm⁻¹). Error denotes the percentage deviation of the retrieved results from the nominal values.

View Table

Based on the measurements at ρ₂ = 60 mm as an example, nominal and retrieved quantities of μ_a and $\mu_{\mathrm{s}}^{\prime}$ of the phantom’s second layer together with their relative errors are summarized in Table 1 for all ten cases where μ_a2 increases gradually with a step of about 0.0185 cm⁻¹. The procedures for determining the nominal values and their uncertainties were presented in Section 3.2. More detailed interpretations are subsequently expatiated in subsection 4.1 and 4.2. The more general analysis of results for the measurements at all ρ₂ are given in subsection 4.3.

4.1 Measured data and results of the first layer

As an example of raw data and the inverse analysis procedures, Fig. 5 shows some measured DTOFs and their corresponding fitting model counterparts from MC simulations. In Fig. 5(a), DTOFs measured at ρ₁ = 10 mm from all 10 cases and the IRF are shown after background subtraction. The DTOFs indicate that an appreciable influence from the second layer on DTOFs only appears in very late time ranges; thus retrieving first layer’s properties by fitting the DTOFs at ρ₁ with the MC homogenous forward model is feasible. Figure 5(b) shows an example of the measured DTOF of case 5 at ρ₁ = 10 mm and its optimal MC homogenous fit (Fit was performed from 50% of the peak at left, to 0.1% at right of the DTOFs. Time 0 is defined as onset of the IRF), and the flat fitting residuals represent high fit quality. By analyzing DTOFs at ρ₁, the optical properties of the first layer are retrieved as (μ_a1, $\mu_{\mathrm{s} 1}^{\prime}$) = (0.100, 9.71) cm⁻¹, with the deviations (−2.0%, +4.0%) from the nominal values (0.102, 9.34) cm⁻¹.

 

Fig. 5. (a) DTOFs measured at ρ₁ = 10 mm for the 10 cases. (b) Fit for case 5: DTOFs at ρ₁ = 10 mm, the optimal fit from homogenous MC and the weighted residuals. (c) DTOFs measured at ρ₂ = 60 mm for the 10 cases. (d) Fit for case 5: DTOFs at ρ₁ = 60 mm, the optimal fit with layered MC and the weighted residuals. (e) DTOFs measured at ρ₂ = 25, 30, 35, 40, and 60 mm (after attenuation filter correction) for case 5, to calculate the SRAC. (f) SRAC for case 5 in selected ρ₂ range, the optimal fit with layered MC and the weighted residuals.

Download Full Size | PDF

Contrarily, as depicted in Fig. 5(c), the prominent effects from the second layer on DTOFs measured at ρ₂ = 60 mm, show that the changes of μ_a2 alter DTOFs throughout their profiles from early time ranges. As one example for the TD inverse process using the MC layered forward model to retrieve optical properties of the second layer, Fig. 5(d) shows the measured and the simulated DTOF for case 5 at ρ₂ = 60 mm. The flat residuals vs. time indicate a very good agreement between both curves. A small (for the case of Fig. 5(d) the minimal) χ² value is obtained, representing the optimal TD solution in {μ_a, $\mu_{\mathrm{s} }^{\prime}$} space. In the analysis here, the TD fits are all performed from 80% of the peak at left, to 0.1% at right of the DTOFs, i.e., almost the whole DTOFs; so all time-resolved information from both layers is considered.

In this work the integrals of photon counts over all time channels, after background subtraction, are used as the SD light intensity information (SRACs). Regions after pulses could be included here thanks to the clean tails of the DTOFs. Figure 5(e) shows DTOFs at all ρ₂ (25, 30, 35, 40, and 60 mm) for case 5. All DTOFs at first 4 ρ₂ are measured with same filter setting while the DTOF at 60 mm is rescaled based on the ratio of filter settings. Note that it is not mandatory to measure all DTOFs with the same filter settings. In practice, appropriate filter ratios and corrections can improve the dynamic range for photon counts. The integrals of DTOFs at all ρ₂ formed the SRAC in Fig. 5(f). After applying the SD inverse model, the optimal layered MC fit can be obtained when the χ² value is minimal.

4.2 Results of the second layer based on different methods

Table 1 provided the overview of the final quantitative results by all three methods, but the comparison of the methods and the trend of the results deserve more attention. In Fig. 6 we present the complete knowledge of the detailed χ² error norm distributions in the {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space under TD, SD, and SeTD method for representative cases 1, 3, 5, 6, 8, and 10. Such error norm surfaces, containing information from TD or SD observations of photon propagation in diffuse media, can be regarded as the temporal or spatial projection of the optical absorption and scattering features on the {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space. To illustrate the relative deviations despite the large span of μ_a2 for the different cases, the horizontal axis (μ_a2) is always given in relative units with respect to the nominal values, while the vertical axis ($\mu_{\mathrm{s} 2}^{\prime}$) is accordingly given in absolute units such that it has the same number of pixels as the horizontal axis. The ratio of $\mu_{\mathrm{s} 2}^{\prime}$ range and μ_a2 range keeps the same for all cases so that details in all the subplots are visible and comparable.

 

Fig. 6. χ² distributions in the {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space. From left to right: TD (ρ = 60 mm), SD, and SeTD; from top to bottom: case 1, 3, 5, 6, 8, and 10. Note that in cases 1 and 3 μ_a1 > μ_a2, in case 5 μ_a1 ≈ μ_a2 and in cases 6, 8 and 10 μ_a1 < μ_a2 . White × crosses are the nominal values, yellow, green, and cyan + crosses are the optimal solutions for TD, SD and SeTD, respectively. Yellow contour lines surround the areas where χ² values are smaller than 2 times minimal χ².

Download Full Size | PDF

Like the results of the simulations in Fig. 1, it can be clearly noticed that the effects from (μ_a2, $\mu_{\mathrm{s} 2}^{\prime}$) properties of the measured phantoms are complementary and symmetric with respect to the χ² error norm surface on TD and SD. As for TD shown in the first column of Fig. 6, μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ generally exhibit a positive correlation, whereas for SD shown in the second column, they are negatively correlated. Such opposite correlations lead to different profiles and orientations of the low-χ² region in the {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space. Although one can still get an optimal solution by both TD and SD methods, i.e. a χ² minimum at one (μ_a2, $\mu_{\mathrm{s} 2}^{\prime}$) combination marked as yellow and green crosses in Fig. 6, they are always quite unstable and sensitive to noise. A small perturbation on measurements may lead to a big deviation from the true values, as implied by the elongated dark red regions (enclosed by yellow contour lines) on the χ² error norm surface of TD and SD.

In the third column of Fig. 6, the above-mentioned complementarity is utilized to carry out the SeTD method by combining TD and SD information, Thereby, more well-defined χ² minima, as the consequence of the better convergence of χ² error norm surface, can always be obtained by the SeTD method to get optimal solutions of the inverse problem. As explained before, effective mutual complementation is crucial to achieving the convergent results.

Based on Fig. 6 and Table 1, several important points should be noted:

4.3 Performance metrics

In order to compare the performance of the SeTD method with respect to the TD and SD methods more quantitatively, two metrics are used to characterize the retrieval of optical properties: error and linearity.

Error: The retrieved values and relative errors, i.e. (retrieved μ / nominal μ) – 1, were shown in Table 1. However, unlike the SD method which always gains the information from all ρ₂, the TD results shown are obtained from the measurements at ρ₂ = 60 mm, and the SeTD results shown include the TD information from ρ₂ = 60 mm as well. According to our analysis, although there is no systematic dependence of the results on ρ₂, the TD measurements at different ρ₂ would lead to varying results for μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ . In other words, the retrieved optical properties from the TD method depend on the location where the measurements took place and thus are unreliable. We considered the retrieved μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ values from measurements at ρ₂ = 25, 30, 35, 40, 60 mm for all 10 cases respectively, and then calculated distributions of the TD and SeTD results at those ρ₂ to draw box plots for μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ and their corresponding relative error with regard to nominal values.

The distributions of retrieved values of μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ from the TD and SeTD methods at different ρ₂ are shown in Figs. 7(a) and 7(b). Clearly distributions from TD are much broader and hence the variances are bigger. Along with the increase of μ_a2, the deviation becomes more apparent. In Figs. 7(c) and 7(d) the relative errors (bias with respect to the nominal values) are presented for the comprehensive investigation between cases of different μ_a2. For μ_a2, the SeTD relative errors clearly stand out for cases 4 to 10. Only for low absorption (cases 1 to 3), where small absolute deviations lead to bigger relative errors, relatively broad distributions are visible in SeTD relative errors. However, low absorption in the deep layer is neither realistic nor of practical interest under in vivo conditions [38]. On the other hand, for $\mu_{\mathrm{s} 2}^{\prime}$ , the results of all cases show that the SeTD can provide more accurate and reliable values than TD at all source-detector distances. In general, the distributions of results from the SeTD method at different ρ₂ are much narrower and the medians of the relative errors are closer to zero, i.e., the SeTD results are largely independent of the source-detector distance at which the TD measurements taken, and have both low variance and low bias with respect to the nominal values.

 

Fig. 7. TD and SeTD comparison of distributions of the results from measurements at different ρ₂ for (a) μ_a2 and (b) $\mu_{\mathrm{s} 2}^{\prime}$ . The corresponding relative error distributions for (c) μ_a2 and (d) $\mu_{\mathrm{s} 2}^{\prime}$ . (e) Linearity of the TD, SD, and SeTD methods and R² distribution over different ρ₂ for TD and SeTD methods. Red lines in the boxes indicate the median, bottom and top edges indicate the 25th and 75th percentiles, and whiskers indicate the minimal and maximal values.

Download Full Size | PDF

Another relevant component of uncertainty of the retrieved values for μ_a2 and $\mu_{\mathrm{s} 2}^{\prime}$ is the uncertainty of the determined parameters of the first layer, i.e. μ_a1, $\mu_{\mathrm{s} 1}^{\prime}$ and its thickness. In order to assess this influence we varied the first layer parameters μ_a1 and $\mu_{\mathrm{s} 1}^{\prime}$ retrieved by the TD analysis at ρ₁ = 10 mm sequentially by ± 9%, and ± 5%, respectively. In two further investigations we considered deviations of the first layer thickness by ± 10%. The Monte Carlo database was extended accordingly. The analyses of the six deviation situations showed that the relative errors of μ_a2 values from SeTD generally remain within ± 10%, with very few exceptions mostly in the cases of lowest and highest μ_a2. Comparatively speaking, SeTD performed clearly better than SD and TD in general. SeTD results surpass those of TD and SD methods for most μ_a2 cases. The μ_a2 relative errors from SD are generally the largest, and those of TD are much less stable. For the situation of + 10% thickness, a pattern with negative relative errors for small μ_a2 cases and positive relative errors for large μ_a2 cases was observed, with a nearly linear trend in between. The pattern was inversed for the situation of - 10% thickness.

For $\mu_{\mathrm{s} 2}^{\prime}$ , the relative errors for all six situations introduced were enlarged for all methods (SD, TD, SeTD). However, the errors from SeTD are smaller than those from TD and SD and are at the same level as the values shown in Fig. 7(d). Having in mind that the most relevant parameter in physiological context is μ_a2, the SeTD method is sufficiently stable with respect to errors in the first layer parameters.

Linearity: Linearity here stands for the linear dependence of retrieved μ against the nominal μ. This metric could provide an effective way to evaluate how well a method can detect the changes in the deep layer without distortions. Also, even if one method retrieves absolute values of μ with an insufficient accuracy, the values of relative changes could still be assessed when the method’s linearity is satisfactory.

We performed a linear fit of the retrieved μ_a2 values versus their nominal values for TD, SD, and SeTD methods, and obtained the linearity by the fit slope and the coefficient of determination R² for the goodness of fit. The linear fits are plotted in Fig. 7(e), in which the μ_a2 results from various ρ₂ are marked as symbols (note that for SD method the symbols represent the results from all ρ₂). The lines are the linear regression for the median of these results from various ρ₂, taking into account uncertainties of μ_a2 derived from the difference of 25th and 75th percentiles for TD and SeTD. Firstly, the slopes of lines from TD, SD, and SeTD methods are 1.02 ± 0.01, 0.89 ± 0.02, and 0.95 ± 0.02, respectively, indicating that in the sense of linearity the SeTD method shows the combination of TD and SD information as well and can comparably detect and maintain the information about relative changes without distortions over a long dynamic range of μ_a2. On the other hand, the median R² values, as shown as the red lines of boxplot in Fig. 7(e), are R² = 0.9872 for TD method and R² = 0.9974 for the SeTD method, showing the more determined linearity of the SeTD method and its ρ₂-independence.

Overall, compared to single-domain methods, the SeTD method can always provide the best performance on error and linearity for retrieval of absorption in the deep layer, and can be independent of the source-detector distance of measurements to a substantial degree. In other words, owing to the better fault tolerance of many conditions during practical measurements, e.g., noise, precise source-detector distance, and accurate superficial layer’s optical properties, the SeTD method can provide more certain and reliable absolute quantity and relative change for deep layer absorption’s assessment.

5. Discussion

The effects of absorption and scattering are entangled in different ways for the measurands in TD and SD. The methodology of the SeTD approach on solving the entanglement, is to integrate the observations from different domains by considering the intrinsic independence nature of absorption and scattering, so as to give the most optimal solution of the ill-conditioned inverse problem.

The SeTD method was experimentally validated by using a two-layered liquid phantom with constant optical properties in the first layer and ten different values of optical properties (primarily absorption) in the second layer. Measurements were carried out at one short source-detector distance and five long distances to distinguish the influences from the first and second layer on diffuse reflectance. The analyses were carried out by TD, SD and the SeTD methods, with Monte-Carlo simulations as the forward model. The retrieved values and error norm distributions on {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space showed that:

If polarization is not considered, in diffuse optics the space-resolved distributions of times of flight of photons represent the maximal information that can be gained from the reflectance scattered field at a single optical wavelength, but also sufficient information for deducing the knowledge of an object’s interior optical absorption and scattering properties when its layered structure is known. In the present work, by taking advantage of the fact that the error norm distributions of TD and SD information are complementary, we combined this information in an artificial spatio-temporal χ_ST² on {μ_a, $\mu_{\mathrm{s}}^{\prime}$} space, to give the optimal quantitative estimate of optical properties in both layers. Instead of measurements requiring a lot of time-domain measurement channels covering a large dynamic range and at various source-detector distances, the SeTD method would only require a short and a long-distance TD measurement combined with a multi-distance SD measurement. It is also worth to note that although in this study DTOFs were recorded at all ρ, it is not a necessary condition for the feasibility of the method. In practice, time-resolved measurements are necessary only at one short and one long distances, while space-resolved amplitude information for the intermediate distances can be accurately obtained by any CW multi-distance detection system. Moreover, rather than isolated fittings between measurements and forward models as single-domain methods focus on, the SeTD method concerns more about the complementary correlations in different domains. The SeTD method takes advantages from and restrains the disadvantages among the domains, to ensure the accuracy and uniqueness of the inverse solution.

In short, for the measured cases, the SeTD method has presented its ability to roughly restrain the retrieval error within 5% for μ_a2 and 10% for $\mu_{\mathrm{s} 2}^{\prime}$ under the condition that the ratio μ_a2/μ_a1 ranges from 0.3 to 2. In previous studies, other researchers have reached similar levels of accuracy on μ_a2 with different approaches, e.g. TD multi-distance measurements with MC simulations [22], TD single distance measurements with optimal estimation algorithm [23], FD multi-distance with diffusion forward model [24], and two-distance TD and SD measurements with an improved diffusion model [25]. In these studies, $\mu_{\mathrm{s} 2}^{\prime}$ was always estimated with large uncertainty due to the aforementioned “$\mu_{\mathrm{s} 2}^{\prime}$ neutrality” phenomenon. To the best of our knowledge, the accuracy level of $\mu_{\mathrm{s} 2}^{\prime}$ achieved here sets a cutting-edge benchmark. Moreover, the SeTD method overcomes some practical problems of single-domain methods, for instance the high uncertainty when measuring at different source-detector distances with the TD method, and the high cross-talk between unknown absorption and scattering with the SD method. The retrieved values by SeTD are independent of the location of the TD detection, and exhibit strong resistance against the cross-talk and the influence of noise owing to the convergence of the error norm. Linearity of SeTD also surpasses single-domain methods, which is crucial for representing changes even if the absolute estimates of optical properties are insufficiently accurate. The SeTD method was also found to be less sensitive for errors in the first layer parameters than the TD and SD methods.

If DTOFs at all relevant ρ₂ have been recorded, a TD global fit could be performed as an alternative to the SeTD method. Such fit contains the SD information in a different manner than SeTD, since it takes the positions of the different DTOFs on the time axis into account. In this way, systematic deviations between experimental and theoretical maxima of the DTOFs could occur. Furthermore, the IRF plays a stronger role than in the SeTD method. The accurate measurement of the IRF is challenging, for it should reflect the illumination of the detection fiber exactly in the same manner as with diffusely reflecting phantoms. The SD information utilized as part of the SeTD method does not depend on the IRF. Hence, the SeTD method overcomes some difficulties related to TD global fit.

There are also limitations of the present study. Firstly, the physical structure of the phantom investigated here is still simplified, although closer to reality than homogeneous ones. It consists of a slab with only 2 layers, and the thickness of the first layer is assumed to be known. Secondly, the optical properties of the first layer have not been changed during the measurements, since we mainly focus on the retrieval performance of the optical properties of the second layer. Finally, despite the absorption in the second layer is investigated for 10 different cases, the scattering is assumed as being homogenous for both layers during the measurements although the values of scattering are considered as unknown. These simplifications may limit the ability of SeTD method to reflect the reality. Further studies will require investigations on complex structures with heterogeneous optical properties.

6. Conclusion

In this study we have developed a new multi-domain approach, namely space-enhanced time domain (SeTD) method, to simultaneously retrieve the absolute quantities of optical absorption and scattering properties in two-layered turbid media. The new methodology of the SeTD method integrates the TD and SD information to give the most optimal estimate of the ill-conditioned inverse problem, rather than isolated fittings of various data types from individual domains. The “deep scattering neutrality” is investigated and summarized in the form of error norm surface. The accuracy, reliability and linearity of the SeTD method are demonstrated to be better than those of the single-domain methods. While the present study was carried out in the near- infrared spectrum, similar methodology and principle could be applied in other spectral ranges so that the method may be used for other general scenarios besides tissues. For the case of measuring at low photon counts, the method also has the potential to obtain more reliable results under the condition of substantial noise. Besides, the concept of merging data from different domains may indicate a new way of solving inverse problems in diffuse optics, namely the integration of different kinds of data types and measurands, such as temporal and central moments of DTOFs, time windows, and frequency domain quantities.