## Abstract

A novel methodology for solving the inverse problem of diffuse optics for two-layered structures is proposed to retrieve the absolute quantities of optical absorption and reduced scattering coefficients of the layers simultaneously. A liquid phantom with various optical absorption properties in the deep layer is prepared and experimentally investigated using the space-enhanced time-domain method. Monte-Carlo simulations are applied to analyze the different measurements in time domain, space domain, and by the new methodology. The deviations of retrieved values from nominal values of both layers’ optical properties are simultaneously reduced to a very low extent compared to the single-domain methods. The reliability and uncertainty of the retrieval performance are also considerably improved by the new methodology. It is observed in time-domain analyses that for the deep layer the retrieval of absorption coefficient is almost not affected by the scattering properties and this kind of “deep scattering neutrality” is investigated and overcome as well.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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^{−1}, while for the 2nd layer the optical properties cover the ranges of *μ*_{a2} ∈ [0.05, 0.25] cm^{−1}, $\mu_{\mathrm{s} 2}^{\prime}$ ∈ [6, 14] cm^{−1}.

The MC simulations of a combination (*μ*_{a2}, $\mu_{\mathrm{s} 2}^{\prime}$) = (0.15, 10) cm^{−1} 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 *χ*^{2} form as the universal objective function:

*m*and

_{i}*s*denote the measurement and the simulation, respectively, with

_{i}*i*= (1…

*n*) representing the effective time channels for TD or

*ρ*channels for SD, of the fitting ranges.

*σ*

_{i}is the standard deviation of

*m*

_{i}. The denominator (

*n*-2) denotes the degrees of freedom of the inverse model applied, given the two fit parameters

*μ*

_{a}, $\mu_{\mathrm{s}}^{\prime}$. For TD as shown in Fig. 1(a), the noise-added DTOF at

*ρ*= 25 mm was selected as measurement, and Eq. (1) was used to compare it with the simulated DTOFs for the entire set of (

*μ*

_{a2}, $\mu_{\mathrm{s} 2}^{\prime}$) combinations in the MC database by calculating

*χ*

_{T}^{2}. For SD as shown in Fig. 1(b), (a) SRAC was prepared in the range

*ρ*∈ [20, 25] mm with step size 1 mm by calculating the integrals of the corresponding noise-added DTOFs, and then compared with its counterparts of all (

*μ*

_{a2}, $\mu_{\mathrm{s} 2}^{\prime}$) combinations in the MC database by calculating

*χ*

_{S}^{2}as well. It is evident that in Figs. 1(a) and 1(b) the

*χ*

^{2}distributions exhibit different orientations, i.e. the low-

*χ*

^{2}regions surrounded by the red contour lines lean towards different directions. These two figures indicate that the opposite correlations of absorption and scattering in TD and SD known from homogeneous cases also exist in the two-layer situation, despite being less prominent. The optimal solutions, where the

*χ*

^{2}values reach the minimum, under the given level of noise, are plotted as green (0.1395, 6.2) cm

^{−1}and yellow (0.1435, 12.04) cm

^{−1}crosses for TD and SD fits and could not be improved by repeating, while the true value (0.15, 10) cm

^{−1}is plotted as white cross.

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}*

^{2}, the ST error norm distribution in the {

*μ*

_{a}, $\mu_{\mathrm{s}}^{\prime}$} space, is introduced as:

*λ*is a Generalized Lagrange multiplier to be optimized by iteration.

*χ*

_{ST}^{2}combines all characteristics of second layer’s optical properties from TD and SD to the ST domain.

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}^{2}* 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

^{−1}, 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.

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^{6}/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*_{0}, 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^{−1}, whereas the second layer’s *μ*_{a2} is varied and investigated in 10 cases in the range [0.03 0.2] cm^{−1} while $\mu_{\mathrm{s} 2}^{\prime}$ is approximately constant at 9.2∼9.3 cm^{−1}. 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.

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*_{1} = 0.05 mm, *S*_{2} = 11.05 mm, *S*_{3} = 0.05 mm, and *S*_{4} = 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^{−1}. 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 *ρ*_{1} = 10 mm as described below. The range of *μ*_{a2} was [0.02 0.22] cm^{−1} with a step size of 0.005 cm^{−1}, and the range of $\mu_{\mathrm{s} 2}^{\prime}$ was [4.5 14.5] cm^{−1} with a step size of 0.25 cm^{−1}. 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 *ρ*_{1} = 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:

- (1) Measuring DTOFs at
*ρ*_{1}= 10 mm; - (2) Measuring DTOFs at
*ρ*_{2}= 25, 30, 35, 40 and 60 mm and obtaining SRACs in this range; - (3) Fitting measured DTOFs at
*ρ*_{1}with the white MC simulations of homogenous forward model (5 cm-thick slab), to retrieve the optical properties of the first layer; - (4) Comparing the MC simulations of layered forward model with the measured DTOFs at a selected
*ρ*_{2}to calculate the*χ*^{2}objective function and obtain the*χ*_{T}^{2}distributions for TD, and with the measured SRACs to calculate and obtain the*χ*_{S}^{2}distributions for SD; - (5) Applying the SeTD method on
*χ*_{ST}^{2}to obtain the optimal (*μ*_{a}, $\mu_{\mathrm{s}}^{\prime}$) for the second layer.

The inverse models for TD and SD were applied by using the *χ*^{2} 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 *χ*^{2} 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 (*ρ*_{1}) and long (*ρ*_{2}) 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^{−1}, $\mu_{\mathrm{s}}^{\prime}$ = 1 cm^{−1}, 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 *ρ*_{1} = 10 mm have travelled in the superficial layer, whereas for *ρ*_{2} = 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.

## 4. Results

Based on the measurements at *ρ _{2}* = 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

^{−1}. 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 ρ

_{2}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 *ρ*_{1} = 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 *ρ*_{1} with the MC homogenous forward model is feasible. Figure 5(b) shows an example of the measured DTOF of case 5 at *ρ*_{1} = 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 *ρ*_{1}, the optical properties of the first layer are retrieved as (*μ*_{a1}, $\mu_{\mathrm{s} 1}^{\prime}$) = (0.100, 9.71) cm^{−1}, with the deviations (−2.0%, +4.0%) from the nominal values (0.102, 9.34) cm^{−1}.

Contrarily, as depicted in Fig. 5(c), the prominent effects from the second layer on DTOFs measured at *ρ*_{2} = 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 *ρ*_{2} = 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) *χ*^{2} 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 *ρ*_{2} (25, 30, 35, 40, and 60 mm) for case 5. All DTOFs at first 4 *ρ*_{2} 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 *ρ*_{2} formed the SRAC in Fig. 5(f). After applying the SD inverse model, the optimal layered MC fit can be obtained when the *χ*^{2} 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 *χ*^{2} 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.

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 *χ*^{2} 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-*χ*^{2} 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 *χ*^{2} 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 *χ*^{2} 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 *χ*^{2} minima, as the consequence of the better convergence of *χ*^{2} 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:

- (1) The yellow contour lines in the 1st column of Fig. 6 enclose the dark red region where
*χ*^{2}values are smaller than 2 times the minimal*χ*^{2}on the entire error norm surface of the TD method. Although the nominal values (white X crosses) always fall into these regions, the retrieved results (yellow crosses) appear in distance from them due to the error norm’s resemblance within low*χ*^{2}region during the inverse process. The flatness within the low*χ*^{2}region indicates that these (*μ*_{a2}, $\mu_{\mathrm{s} 2}^{\prime}$) combinations in the region are practically identical for the solver and implies the non-uniqueness of the TD solutions. - (2) The SD
*χ*^{2}error norm surface in the 2nd column of Fig. 6 is comparably much steeper. And consequently, the yellow contour lines enclose a smaller and more confined region with convergent low*χ*^{2}. However, the nominal values hardly fall into these regions, which makes the simultaneous recovery of absorption and scattering nearly impossible. The retrieved results (green crosses) still appear in some distance from nominal values. At the cost of the underestimating $\mu_{\mathrm{s} 2}^{\prime}$, the retrieval performance of*μ*_{a2}is quite unstable among the cases, which is shown in the third column of Table 1 as well. - (3) The TD low
*χ*^{2}regions, i.e. the red flat strips of the error norm surface in Fig. 6, are roughly parallel to the $\mu_{\mathrm{s} 2}^{\prime}$ axis, which denotes that $\mu_{\mathrm{s} 2}^{\prime}$ has very little influence on the DTOF shape of the diffuse reflectance. Namely, it is nearly impossible to determine $\mu_{\mathrm{s} 2}^{\prime}$ through TD information alone. Such uncertainty also manifests itself as the big deviation of the TD retrieved $\mu_{\mathrm{s} 2}^{\prime}$ from nominal values shown in the second column of Table 1. For the cases which*μ*_{a2}values are not far from*μ*_{a1}values, it is particularly apparent that even if the retrieved $\mu_{\mathrm{s} 2}$ is inaccurate, it is nonetheless still possible to obtain*μ*_{a2}with reasonable accuracy by the TD method. While for the cases with stronger contrast between layers, this “$\mu_{\mathrm{s} 2}^{\prime}$ neutrality” declines to a certain extent and the ambiguity of*μ*_{a2}and $\mu_{\mathrm{s} 2}^{\prime}$ increases again. - (4) The SD low
*χ*^{2}regions, i.e. the elongated “valley” enclosed by yellow contour lines, are turning towards the vertical direction with increasing*μ*_{a2}. The trend of retrieval error of*μ*_{a2}in Table 1 indicates that the SD method is more sensitive to high absorption in the second layer. However, it is important to note here that the*χ*^{2}minimum cannot be well allocated due to the occurrence of several local minima along the valley center line caused by noise. Furthermore, the results for the SD method in Table 1 and Fig. 6 are based on the assumption of Poisson noise. Repeated measurements show that the noise of the SD intensities is larger and amounts to about 1%, most likely due to fluctuations of the laser intensity. Taking the increased noise into account we observe a shift of the*χ*^{2}minimum along the valley center line towards the upper left direction. Overall, the stable information from the SD method is that all combinations of*µ*_{a}and $\mu_{\mathrm{s}}^{\prime}$ along the valley center line are consistent with the experimental data. This line can be interpreted, at least qualitatively, as an effective attenuation, defined as*μ*_{eff}= (3*μ*_{a}$\mu_{\mathrm{s}}^{\prime}$)^{1/2}. Generally, the non-uniqueness of the SD method discussed in [17] takes effect here as well. - (5) Likewise to the SD method, along with increasing absorption in the second layer, the low-
*χ*^{2}regions become broader for TD method, too. This is due to the relatively fewer collected photons in the later time ranges, which have experienced the second layer and then survived to the exit surface. In general, lower signal-to-noise ratio (SNR) will always cause a higher uncertainty. Under lower SNR conditions (low photon counts or high noise), the retrieval accuracy for both absorption and scattering would be reduced, since an error in one unknown would manifest itself as an artefact to the other. - (6) As can be seen from the bold values in Table 1, SeTD reached the best performance among the 3 methods for almost all cases. It always derives the $\mu_{\mathrm{s} 2}^{\prime}$ values closest to the nominal values. The accuracy is also improved for
*μ*_{a2}when comparing with original results from the TD method. A well-defined convergence in the {*μ*_{a}, $\mu_{\mathrm{s}}^{\prime}$} space is obtained for all cases by SeTD. As illustrated in the 3rd column of Fig. 6, the SeTD method maintains the “$\mu_{\mathrm{s} 2}^{\prime}$ neutrality” from TD and, at the same time, benefits of the steepness and convergence in the low*χ*^{2}regions from SD. With the SeTD method, the advantages from both TD and SD are effectively combined through mutual complementation.

*χ*

^{2}region can be reduced, but also the crosstalk among

*μ*

_{a2}and $\mu_{\mathrm{s} 2}^{\prime}$ can be restrained. The entangled contributions from scattering and absorption are decoupled in this artificial spatio-temporal projection

*χ*

_{ST}^{2}by integrating SD and TD information.

#### 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 *ρ*_{2}, the TD results shown are obtained from the measurements at *ρ*_{2} = 60 mm, and the SeTD results shown include the TD information from *ρ*_{2} = 60 mm as well. According to our analysis, although there is no systematic dependence of the results on *ρ*_{2}, the TD measurements at different *ρ*_{2} 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 *ρ*_{2} = 25, 30, 35, 40, 60 mm for all 10 cases respectively, and then calculated distributions of the TD and SeTD results at those *ρ*_{2} 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 *ρ*_{2} 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 *ρ*_{2} 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.

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 *ρ*_{1} = 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*^{2} for the goodness of fit. The linear fits are plotted in Fig. 7(e), in which the *μ*_{a2} results from various *ρ*_{2} are marked as symbols (note that for SD method the symbols represent the results from all *ρ*_{2}). The lines are the linear regression for the median of these results from various *ρ*_{2}, 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*^{2} values, as shown as the red lines of boxplot in Fig. 7(e), are *R*^{2} = 0.9872 for TD method and *R*^{2} = 0.9974 for the SeTD method, showing the more determined linearity of the SeTD method and its *ρ*_{2}-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:

- (1) For TD method, the error norm is sensitive to the change of absorption in the second layer, while it is virtually insensitive to scattering in the second layer, i.e., the temporal profiles of diffuse reflectance are practically insensitive to $\mu_{\mathrm{s} 2}^{\prime}$ . As introduced above by the term of “deep scattering neutrality”, it is nearly impossible to simultaneously recover
*μ*_{a2}and $\mu_{\mathrm{s} 2}^{\prime}$ through TD alone, if no*prior*knowledge is provided. - (2) For SD method, under the presence of complex noise, the retrieval performance is bad and the uncertainty is high. Due to the limited information one can acquire from SD, the solutions of
*μ*_{a2}and $\mu_{\mathrm{s} 2}^{\prime}$ with the most confidence, i.e. the low-*χ*^{2}region on the error norm distribution, always appear at some distance from the nominal values, indicating that it is hard to retrieve absorption and scattering in the second layer through the SD method alone. - (3) For both TD and SD methods, there is no clear convergence on the error norm distribution for the second layer’s properties, indicating that the ambiguity within the low-
*χ*^{2}region of error norm distributions would make the retrieval unstable. Conversely, for the SeTD method, better convergence of the error norm can be achieved, which yields a well-defined and reliable solution. In all measured cases, the*μ*_{a2}accuracy is improved in comparison to the results from conventional TD, and the most accurate $\mu_{\mathrm{s} 2}^{\prime}$ values are always from the SeTD method as well.

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}*

^{2}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 *ρ _{2}* 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.

## Funding

H2020 Marie Skłodowska-Curie Actions (BitMap, Grant No. 675332, Innovative Training Networks); Narodowe Centrum Nauki (UMO-2019/33/N/ST7/02918).

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **Z. A. Steelman, D. S. Ho, K. K. Chu, and A. Wax, “Light-scattering methods for tissue diagnosis,” Optica **6**(4), 479–489 (2019). [CrossRef]

**2. **A. M. Smith, M. C. Mancini, and S. Nie, “Bioimaging: second window for in vivo imaging,” Nat. Nanotechnol. **4**(11), 710–711 (2009). [CrossRef]

**3. **F. F. Jöbsis-VanderVliet, “Discovery of the near-infrared window into the body and the early development of near-infrared spectroscopy,” J. Biomed. Opt. **4**(4), 392–397 (1999). [CrossRef]

**4. **G. Bale, C. E. Elwell, and I. Tachtsidis, “From Jöbsis to the present day: a review of clinical near-infrared spectroscopy measurements of cerebral cytochrome-c-oxidase,” J. Biomed. Opt. **21**(9), 091307 (2016). [CrossRef]

**5. **F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science **198**(4323), 1264–1267 (1977). [CrossRef]

**6. **T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. **73**(7), 076701 (2010). [CrossRef]

**7. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**(2), R41–R93 (1999). [CrossRef]

**8. **S. R. Arridge and J. C. Schotland, “Optical tomography: Forward and inverse problems,” Inverse Probl. **25**(12), 123010 (2009). [CrossRef]

**9. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,” Appl. Opt. **28**(12), 2331–2336 (1989). [CrossRef]

**10. **F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. M. Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” NeuroImage **85**, 6–27 (2014). [CrossRef]

**11. **D. J. Cuccia, F. Bevilacqua, A. J. Durkin, and B. J. Tromberg, “Modulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain,” Opt. Lett. **30**(11), 1354–1356 (2005). [CrossRef]

**12. **D. Lighter, J. Hughes, I. Styles, A. Filer, and H. Dehghani, “Multispectral, non-contact diffuse optical tomography of healthy human finger joints,” Biomed. Opt. Express **9**(4), 1445–1460 (2018). [CrossRef]

**13. **J. B. Fishkin and E. Gratton, “Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge,” J. Opt. Soc. Am. A **10**(1), 127–140 (1993). [CrossRef]

**14. **L. Wang, S. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Computer methods and programs in biomedicine **47**(2), 131–146 (1995). [CrossRef]

**15. **A. Liemert and A. Kienle, “Infinite space Green’s function of the time-dependent radiative transfer equation,” Biomed. Opt. Express **3**(3), 543–551 (2012). [CrossRef]

**16. **A. Liemert, D. Reitzle, and A. Kienle, “Analytical solutions of the radiative transport equation for turbid and fluorescent layered media,” Sci. Rep. **7**(1), 3819 (2017). [CrossRef]

**17. **S. R. Arridge and W. R. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. **23**(11), 882–884 (1998). [CrossRef]

**18. **M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. **44**(7), 1699–1717 (1999). [CrossRef]

**19. **M. Wolf, G. Naulaers, F. van Bel, S. Kleiser, and G. Greisen, “A review of near infrared spectroscopy for term and preterm newborns,” J. Near Infrared Spectrosc. **20**(1), 43–55 (2012). [CrossRef]

**20. **B. Hallacoglu, A. Sassaroli, S. Fantini, and A. Troen, “Cerebral perfusion and oxygenation are impaired by folate deficiency in rat: absolute measurements with noninvasive near-infrared spectroscopy,” J. Cereb. Blood Flow Metab. **31**(6), 1482–1492 (2011). [CrossRef]

**21. **B. Hallacoglu, A. Sassaroli, S. Fantini, M. Wysocki, E. Guerrero-Berroa, M. Beeri, V. Haroutunian, M. Shaul, I. Rosenberg, and A. Troen, “Absolute measurement of cerebral optical coefficients, hemoglobin concentration and oxygen saturation in old and young adults with near-infrared spectroscopy,” J. Biomed. Opt. **17**(8), 081406 (2012). [CrossRef]

**22. **J. Selb, T. M. Ogden, J. Dubb, Q. Fang, and D. A. Boas, “Comparison of a layered slab and an atlas head model for Monte Carlo fitting of time-domain near-infrared spectroscopy data of the adult head,” J. Biomed. Opt. **19**(1), 016010 (2014). [CrossRef]

**23. **F. Martelli, S. Del Bianco, L. Spinelli, S. Cavalieri, P. Di Ninni, T. Binzoni, A. Jelzow, R. Macdonald, and H. Wabnitz, “Optimal estimation reconstruction of the optical properties of a two-layered tissue phantom from time-resolved single-distance measurements,” J. Biomed. Opt. **20**(11), 115001 (2015). [CrossRef]

**24. **B. Hallacoglu, A. Sassaroli, and S. Fantini, “Optical characterization of two-layered turbid media for non-invasive, absolute oximetry in cerebral and extracerebral tissue,” PLoS One **8**(5), e64095 (2013). [CrossRef]

**25. **A. Kienle, T. Glanzmann, G. Wagnieres, and H. van den Bergh, “Investigation of two-layered turbid media with time-resolved reflectance,” Appl. Opt. **37**(28), 6852–6862 (1998). [CrossRef]

**26. **M. Shimada, C. Sato, Y. Hoshi, and Y. Yamada, “Estimation of the absorption coefficients of two-layered media by a simple method using spatially and time-resolved reflectances,” Phys. Med. Biol. **54**(16), 5057–5071 (2009). [CrossRef]

**27. **L. Yang, H. Wabnitz, T. Gladytz, R. Macdonald, and D. Grosenick, “Spatially-enhanced time-domain NIRS for accurate determination of tissue optical properties,” Opt. Express **27**(19), 26415–26431 (2019). [CrossRef]

**28. **S. Suzuki, S. Takasaki, T. Ozaki, and Y. Kobayashi, “Tissue oxygenation monitor using NIR spatially resolved spectroscopy,” Proc. SPIE **3597**, 582–592 (1999). [CrossRef]

**29. **Q. Fang and D. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express **17**(22), 20178–20190 (2009). [CrossRef]

**30. **A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A **29**(10), 2110–2117 (2012). [CrossRef]

**31. **C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. **18**(5), 050902 (2013). [CrossRef]

**32. **E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. **13**(6), 060504 (2008). [CrossRef]

**33. **P. Di Ninni, F. Martelli, and G. Zaccanti, “Intralipid: towards a diffusive reference standard for optical tissue phantoms,” Phys. Med. Biol. **56**(2), N21–N28 (2011). [CrossRef]

**34. **P. Di Ninni, F. Martelli, and G. Zaccanti, “The use of India ink in tissue simulating phantoms,” Opt. Express **18**(26), 26854–26865 (2010). [CrossRef]

**35. **L. Spinelli, M. Botwicz, N. Zolek, M. Kacprzak, D. Milej, P. Sawosz, A. Liebert, U. Weigel, T. Durduran, F. Foschum, A. Kienle, F. Baribeau, S. Leclair, J.-P. Bouchard, I. Noiseux, P. Gallant, O. Mermut, A. Farina, A. Pifferi, A. Torricelli, R. Cubeddu, H.-C. Ho, M. Mazurenka, H. Wabnitz, K. Klauenberg, O. Bodnar, C. Elster, M. Bénazech-Lavoué, Y. Bérubé-Lauzière, F. Lesage, D. Khoptyar, A. A. Subash, S. Andersson-Engels, P. Di Ninni, F. Martelli, and G. Zaccanti, “Determination of reference values for optical properties of liquid phantoms based on Intralipid and India ink,” Biomed. Opt. Express **5**(7), 2037–2053 (2014). [CrossRef]

**36. **F. Martelli, P. Di Ninni, G. Zaccanti, D. Contini, L. Spinelli, A. Torricelli, R. Cubeddu, H. Wabnitz, M. Mazurenka, R. Macdonald, A. Sassaroli, and A. Pifferi, “Phantoms for diffuse optical imaging based on totally absorbing objects, part 2: experimental implementation,” J. Biomed. Opt. **19**(7), 076011 (2014). [CrossRef]

**37. **L. Spinelli and G. Zaccanti, “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. time-resolved method,” Opt. Express **15**(11), 6589–6604 (2007). [CrossRef]

**38. **A. Farina, A. Torricelli, I. Bargigia, L. Spinelli, R. Cubeddu, F. Foschum, M. Jäger, E. Simon, O. Fugger, A. Kienle, F. Martelli, P. Di Ninni, G. Zaccanti, D. Milej, P. Sawosz, M. Kacprzak, A. Liebert, and A. Pifferi, “In-vivo multilaboratory investigation of the optical properties of the human head,” Biomed. Opt. Express **6**(7), 2609–2623 (2015). [CrossRef]

**39. **S. Del Bianco, F. Martelli, F. Cignini, G. Zaccanti, A. Pifferi, A. Torricelli, A. Bassi, P. Taroni, and R. Cubeddu, “Liquid phantom for investigating light propagation through layered diffusive media,” Opt. Express **12**(10), 2102–2111 (2004). [CrossRef]

**40. **E. Alerstam, S. Andersson-Engels, and T. Svensson, “White Monte Carlo for time-resolved photon migration,” J. Biomed. Opt. **13**(4), 041304 (2008). [CrossRef]

**41. **S. Feng, F. Zeng, and B. Chance, “Photon migration in the presence of a single defect: a perturbation analysis,” Appl. Opt. **34**(19), 3826–3837 (1995). [CrossRef]

**42. **A. Liebert, H. Wabnitz, N. Zołek, and R. Macdonald, “Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media,” Opt. Express **16**(17), 13188–13202 (2008). [CrossRef]