Time-resolved near infrared light propagation using frequency domain superposition

Time-resolved temporal point spread function (TPSF) measurement of near infrared spectroscopic (NIRS) data allows the estimation of absorption and reduced scattering properties of biological tissues. Such analysis requires an iterative calculation of the theoretical TPSF curve using mathematical and computational models of the domain being imaged which are computationally complex and expensive. In this work, an efficient methodology for representing the TPSF data using a superposition of cosines calculated in frequency domain is presented. The proposed method is outlined and tested on finite element realistic models of the human neck and head. Using an adult head model containing ~140k nodes, the TPSF calculation at each node for one source is accelerated from 3.11 s to 1.29 s within an error limit of ± 5% related to the time domain calculation method. © 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (170.5270) Photon density waves; (170.3660) Light propagation in tissues. References and links 1. A. Pifferi, D. Contini, A. D. Mora, A. Farina, L. Spinelli, and A. Torricelli, “New frontiers in time-domain diffuse optics, a review,” J. Biomed. Opt. 21(9), 091310 (2016). 2. A. Torricelli, D. Contini, A. Pifferi, M. Caffini, R. Re, L. Zucchelli, and L. Spinelli, “Time domain functional NIRS imaging for human brain mapping,” Neuroimage 85(Pt 1), 28–50 (2014). 3. H. Wabnitz, M. Moeller, A. Liebert, H. Obrig, J. Steinbrink, and R. Macdonald, “Time-resolved near-infrared spectroscopy and imaging of the adult human brain,” Adv. Exp. Med. Biol. 662, 143–148 (2010). 4. W. Weigl, D. Milej, D. Janusek, S. Wojtkiewicz, P. Sawosz, M. Kacprzak, A. Gerega, R. Maniewski, and A. Liebert, “Application of optical methods in the monitoring of traumatic brain injury: A review,” J. Cereb. Blood Flow Metab. 36(11), 1825–1843 (2016). 5. M. Kacprzak, A. Liebert, W. Staszkiewicz, A. Gabrusiewicz, P. Sawosz, G. Madycki, and R. Maniewski, “Application of a time-resolved optical brain imager for monitoring cerebral oxygenation during carotid surgery,” J. Biomed. Opt. 17(1), 016002 (2012). 6. D. Grosenick, H. Rinneberg, R. Cubeddu, and P. Taroni, “Review of optical breast imaging and spectroscopy,” J. Biomed. Opt. 21(9), 091311 (2016). 7. D. A. Boas, L. E. Campbell, and A. G. Yodh, “Scattering and Imaging with Diffusing Temporal Field Correlations,” Phys. Rev. Lett. 75(9), 1855–1858 (1995). 8. T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” Neuroimage 85(Pt 1), 51–63 (2014). 9. C. Lindner, M. Mora, P. Farzam, M. Squarcia, J. Johansson, U. M. Weigel, I. Halperin, F. A. Hanzu, and T. Durduran, “Diffuse Optical Characterization of the Healthy Human Thyroid Tissue and Two Pathological Case Studies,” PLoS One 11(1), e0147851 (2016). 10. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41 (1999). 11. A. Liebert, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, and H. Rinneberg, “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. 42(28), 5785–5792 (2003). 12. A. Liebert, H. Wabnitz, and C. Elster, “Determination of absorption changes from moments of distributions of times of flight of photons: optimization of measurement conditions for a two-layered tissue model,” J. Biomed. Opt. 17(5), 057005 (2012). 13. A. Pifferi, P. Taroni, G. Valentini, and S. Andersson-Engels, “Real-time method for fitting time-resolved reflectance and transmittance measurements with a monte carlo model,” Appl. Opt. 37(13), 2774–2780 (1998). Vol. 9, No. 1 | 1 Jan 2018 | BIOMEDICAL OPTICS EXPRESS 41


Introduction
Time-resolved spectroscopy (TRS) devices [1][2][3] measure a temporal point spread function (TPSF) at several wavelengths within the near infrared window (650 -1000 nm).The TPSF is also called a distribution of time of flight of photons (DTOF), which represents the DTOF measured in tissue at some source and detector distances using a pulsed (pico or femto seconds in width) light source and single photon counting time-resolved detectors placed usually on the tissue surface in reflectance or transmission geometry.Shape of the measured TPSF curve depends on both the often heterogeneous absorption μa and reduced scattering μ′ s properties of the medium.Thus, analysis the TPSF shape allow the estimation of the optical properties at different wavelengths to allow recovery of concentrations of chromophores e.g.oxygenated and reduced hemoglobin.
TRS has been used in many in-vivo and clinical studies.In particular it was used to assess condition of traumatic brain injury patients [4], brain oxygenation during carotid surgery [5] and for optical mammography [6] or functional brain mapping [2].TRS has also been combined with diffuse correlation spectroscopy (DCS) [7,8] for noninvasive recovery of hemodynamic parameters of human thyroid on healthy subjects and cancer patients [9], where the TRS was utilized to measure thyroid oxygenation and the optical properties needed to assess thyroid perfusion with the DCS technique.In all of these cases, the calculation and estimation of the underlying optical properties relies on some model-based parameter fitting algorithm, which can be computationally expensive [10].
To allow the fitting of the optical parameters to the measured data, the shape of the TPSF can be analyzed using different approaches, ranging from extracted curve parameters to full shape analysis.Using the normalized statistical moments of the curve where the number of photons (intensity), mean time of flight and variance can be used to estimate bulk or twohomogeneous layers absorption μa and reduced scattering μ′ s properties of the medium [11,12].Further, an iterative curve fitting procedure can be applied to improve this recovery [13], which minimizes the error between the measured and modelled TPSF with respect to optical properties of the medium.Others have also proposed that Mellin-Laplace transform of the TPSF can also be used for tomographic reconstruction of tissue absorption [14], specifically when concerned with changes in deep tissue where the so called 'late-photons' carry most of the information.
The curve fitting procedure requires calculation of time-resolved data using a computational model of light propagation for given optical properties in an iterative manner, which as a consequence requires a fast and reliable model for a complex heterogeneous system [15,16].
The gold standard of TPSF generation is the Monte Carlo (MC) method [17] as applied in for example layered models [18] or arbitrary heterogeneous structures represented by voxels or finite element meshes [19].Despite the accuracy of the MC (the best physical representation of photon transport within a turbid medium [20]) and latest improvements on speed (graphical processing unit (GPU) accelerated MC) [21] the MC generated TPSFs will always be inherently slow and noisy, particularly when considering late photons, large source/detector separation or highly attenuating material.MC is a stochastic process where the TPSF noise amplitude approaches to zero when the calculation time (number of simulated photons) goes to infinity.Further, fast MC generation of TPSF for a complex model requires finely tuned parallel code executed on state of the art GPUs or a cluster of GPUs.It has been reported that GPU accelerated MC code executed on a single GPU will generate the timeresolved photon fluence rate within MRI based head model (1 mm 3 voxel size) in about 1.7 min for 10 7 photons [21].MC method can be used to calculate a lookup table with TPSF discretized against absorption and scattering.However, the resolution of recovered optical properties will always be limited and the lookup table should be generated for each tissue model separately.Thus, the lookup table approach is mainly utilized for semi-infinite or homogeneously layered tissue models.
Another approach for TPSF generation is an analytical method based on diffusion approximation of photon transport equation [22].The time-resolved diffusion equation benefits from simple analytical solutions for homogenous semi-infinite or infinite medium [23].Alternatively, a complex tissue model can be described with finite element models (FEM) where the diffusion equation can be numerically solved for all degrees of freedom (mesh nodes representing spatial distribution) [24,25] as implemented in NIRFAST software package used in this study (www.nirfast.org).However, this approach for time-resolved curves requires solving a differential equation at many time points representing the TPSF [22,26], typically implemented using the Crank-Nicolson approximation.In this approach, the number of time points define the accuracy of the TPSF and calculation speed.
As an alternative approach to calculating the photon fluence as a function of time, Kienle et al. in [27,28] used a superposition of series of solutions of diffusion equation in frequency domain to generate the time-resolved TPSF.This method was used for analytical solution in two-layered model and to date has not been investigated for more complex heterogeneous models based on numerical solutions.Here, it is observed that this frequency domain approach can be extended for complex finite element models and with appropriately limited frequency band used in the superposition, it may provide one of the computationally fastest approaches in generating time-resolved data within complex heterogeneous models.
In this work, we present the methodology of representing TPSF by superposition of cosine functions calculated in frequency domain to generate the time-resolved TPSF.The method has been tested on heterogeneous finite element models of the human neck and head showing negligible loss of accuracy and a significant gain in calculation speed, as compared to conventional calculation in time domain.
This work is motivated by a need of fast and reliable method of simulation of TPSFs for realistic human neck model and two-layered model in order to recover thyroid optical properties and oxygenation [9].This algorithm has been designed to be a part of analysis suite of the Laser and Ultrasound Co-Analyzer for Thyroid Nodules developed under the European project (LUCA, http://www.luca-project.eu)-the purpose of which is building a multi-modal instrument for thyroid cancer screening.

Methodology
In this section the theory of representation of time-resolved time point spread function by series of simulations in the frequency domain is presented.Additionally, a practical approach of choosing frequencies needed to accurately represent the time-resolved data is shown.

Theory
The time-resolved diffusion approximation of photon transport equation in an optically turbid medium represented by a finite element model can be expressed as: where ) is photon fluence rate (TPSF), c m (r) is speed of light in the medium which is a function of the medium's refractive index n(r), r s is source location and S(r s ,t) is a source term e.g. a pulse of light (Dirac delta δ(r s ,t)) at time t=0.The TPSF Φ(r,t) can be calculated by numerically solving a set of differential equations as is currently implemented in NIRFAST software [22,26]: Equation ( 2) can be numerically solved within a time span t ∈ [0; T] for N t points, where the source term S(r s ,t) represents a pulse at t = 0 and amplitude 1.
The Dirac delta light source can be represented by superposition of infinite number of cosines of the same amplitude A and oscillating at infinite number of angular frequencies ω.For a single frequency ω, Eq. ( 1) can be transformed into the frequency domain using Fourier transform as follows: where the Fourier transform  of the source term δ(r s ,t) is equal to A at frequency ω and 0 elsewhere.Using the superposition principle, Eq. ( 3) can be solved for an infinite number of sources (for infinite number of frequencies) where the right hand side will always equal A.
Finally, infinite number of solutions in the frequency domain Φ(r,ω) can be transformed back into the time domain using inverse Fourier transform as follows: Discrete representation of Eq. ( 4) requires reducing number of frequencies from infinity to some natural number N ω .

Frequencies needed to represent time-resolved data
Contribution of the Φ(r,ω) to the TPSF Φ(r,t) (Eq.( 4)) varies with frequency ω and at some frequencies this contribution may be considered as negligible.Thus, the infinite angular frequency range ω ∈ [0; ∞) may be limited with a minor loss of information carried by photons traveling from the source to the detector.The lowest usable frequency ω 0 can be determined by the maximum observable photons travel time T. For example in time-resolved diffuse optical spectroscopy and for biological samples, this time is measured within a few nanoseconds.In practice, the time scale of TRS hardware [29] for most clinical applications such as human brain and breast imaging (source/detector separation of approximately 7-9 cm [30]) is usually t ∈ [0; 10] ns .The maximum observation time T is always assumed and estimated a-priori for all approaches of TPSF simulations, e.g.MC. Figure 1 shows a schematic representation of contributions to a TPSF curve of angular frequencies close to ω = π/(2T).As shown in Fig. 1a the contributions of frequency components at frequencies ω < π/(2T) approach a constant value A for the entire observation window t ∈ [0; T].Thus, these frequencies of ω < π/(2T) will have a low influence on the TPSF curve shape.However, the TPSF curve is always shifted to the right in the time axis, for source-detector separations of r sd > 0 as there are no instantly detected photons and the time/phase shift for frequency components of low ω will correspond to the maximum of the TPSF curve in such instances.Thus, the maximum of the TPSF curve will be located within the range of t ∈ (0; T) which further weakens the influence of frequencies ω < π/(2T) as shown in Fig. 1b.Considering the above, it is assumed that the lowest usable frequency should be set at ω 0 = π/(2T).
The highest usable frequency ω N can be defined using Nyquist theorem.That is, if the number of frequencies used to reproduce the TPSF curve is defined by a natural number N ω , the maximum usable frequency can be defined as ω N = 2N ω ω 0 = N ω π/T.Finally, Eq. ( 4) can then be approximated at a point in time t by the discrete form as: where ω k = k(ω Nω 0 )/N ω + ω 0 , ω 0 = π/(2T), ω N = 2N ω ω 0 and N ω -number of evenly spaced frequencies within [ω 0 ; ω N ], t ∈ [0; ∞).The unit of fluence in Eq. ( 5) is mm −2 s −1 .The value of N ω can be used to control relation between speed and accuracy of generating TPSFs with Eq. ( 5).This relation will be investigated in the next section.

Validation
In this section the new frequency based calculation of the TPSF (Eq.( 5)) is evaluated against the traditional method (Eq.( 2) -time domain) as already implemented in NIRFAST using a finite element models of homogenous slab, a human head and a human neck.The test was performed on a PC equipped with 10-core Intel Xeon E5-3660 v3 CPU working @ 2.6 GHz and NVidia GeForce GTX 1080 GPU with 2560 CUDA parallel processor cores.Execution times were averaged for 10 runs for each case.
Equation ( 2) and ( 3) can be expressed as system of linear equations AΦ = q as shown in [22,25,26], where A N × N is a sparse matrix representing finite element model defined by N vertex nodes in R 3 Euclidean space, q N × M is composed of M light source vectors and Φ N × M is the unknown photon fluence rate for M sources.In this work, the sparse matrices A and q are transferred to GPU memory once where the problem is solved for each source independently.Thus, the execution time per source decreases exponentially with number of simulated sources M. The execution time per source are shown where all calculations were carried out for M = 8 sources.

Finite element models
The homogenous slab (120 × 120 × 50 mm) is represented by two finite element models of low (LO) and high (HI) resolution meshes quantified as mean volume of tetrahedral FEM elements as shown in Table 1.All meshes used have a resolution that is uniform within the model.Optical properties of both the homogenous slabs LO and HI were set as: μa = 0.01 mm −1 , μ′ s = 1 mm −1 .
The neck model is developed using images registered with 1.5T Siemens magnetic resonance imager using echo planar imaging (T2-weighted) sequence (matrix 256 × 256, resolution 0.84 × 0.84 mm, slice thickness 0.9 mm, 88 slices).The neck volume represented by 88 grey-scale images has been semi-automatically segmented into 8 structures using NIRFASTSlicer [24] software.For the head model, MRI data from a given subject is used together with Statistical Parametric Mapping (SPM) [31,32] which first allows a parametric segmentation of the 5 tissue types (scalp, skull, cerebrospinal fluid, gray and white matter) based on the pixel intensity probability function distribution.The segmented neck and head volumes were transformed into finite element mesh of linear tetrahedrons, carried out with the NIRFASTSlicer meshing module based on Computational Geometry Algorithms Library [35].Optical properties of the human neck and head models at 780 nm are shown in Table 2. Reduced scattering coefficient of cerebrospinal fluid was assessed by inverse of its thickness and set to 0.3 mm −1 as shown in [36].This is about two orders of magnitude higher than values used in Monte Carlo methods [37,38].However, as argued in [36], μ′ s = 0.3 mm −1 allows accurate light propagation modeling and satisfies diffusion approximation condition μ′ s  μ a .The refractive index n=1.33 is homogenous on all models.Quality/resolution of all finite element models is shown in Table 1.
For each model, a grid of 8 sources and 9 detectors are placed at each model surface as shown in Fig. 2. The inter-optode distance within the grid is 13 mm.

Accuracy and speed
Numerical solvers for the time-resolved forward problem as defined in Eq. ( 2) use a solution of the previous time points t n-s (s ≥ 1) to solve the problem at point t n .Thus, the solution accuracy increases for shorter integration steps Δt = t n -t n-1 = T/N t (assuming evenly spaced time samples).The relationship between the accuracy and stability of the solution in time domain with respect to the number of time points N t = 2 k , k ∈ {6; 7; …; 10} is presented in Fig. 3 assuming the neck model (Fig. 2c), observation time of T = 2 ns and source-detector separations 29 mm and 44 mm.The relative residuals are calculated with respect to the most accurate TPSF generated for N t = 2 10 .The accuracy for the neck model is shown as an example since similar findings were found for the head model.The TPSF region of interest (ROI) is defined as a percentage of the maximum value of a TPSF, which is set to 80% on the leading edge and 5% on the falling edge as shown in Fig. 3 as it has been reported [40][41][42][43] that the 50%-100% leading edge limit and 20% -1% falling edge limit can be used to carry out the procedure of curve fitting.Fig. 3. Relation between accuracy and number of time points N t of NIRFAST time domain method as defined in Eq. ( 2).Data shown for the neck model as in Fig. 2c and two sourcedetector pairs.r sd -source-detector separation of 29 and 44 mm.80% and 5% mark TPSF region of interest calculated as a percentage of the maximum value of TPSF for N t = 1024.Residuals are related to TPSFs calculated for N t = 1024.
Accuracy of the TPSF calculated in the frequency domain Eq. ( 5) depends on the number of frequencies N ω (between the lowest and highest outlined above) used to build the curve.The relationship between accuracy and number of frequencies N ω = 2 k , k ∈ {5; 6; …; 9} is presented in Fig. 4 assuming the neck model (Fig. 2c), observation time T = 2 ns, which is always discretized into 1024 samples.The relative residuals are calculated in respect to the most accurate TPSF generated in time domain for N t = 1024.The lowest usable frequency needed to reproduce the TPSF is set to f 0 = ω 0 /(2π) = 1/(4T) which is equal to 125 MHz at T = 2 ns.The upper limit of the frequency band changes with number of frequencies N ω as follows: f N = ω N /(2π) = N ω /(2T) which is equal to {8; 16; 32; 64; 128} GHz at T=2 ns.Fig. 4. Relation between accuracy and number of frequencies N ω of NIRFAST frequency domain method as defined in Eq. ( 5).Data shown for the neck model as in Fig. 2(c) and two source-detector pairs.r sd -source-detector separation of 29 and 44 mm.80% and 5% mark TPSF region of interest calculated as a percentage of the maximum value of TPSF calculated in time domain for N t = 1024 Residuals are related to TPSFs calculated in time domain for N t = 1024.
Accuracy of both methods Eq. ( 2) and ( 5) is compared with analytical solution [39] of diffusion equation assuming semi-infinite medium and Robin boundary condition to match NIRFAST implementation [25].Only top surface of the FEM slab models is considered as a boundary.Figure 5 shows maximum and mean of absolute values of the relative residuals within the ROI for all unique source-detector separations r sd as in Fig. 2a.In case of the slab models the observation time is T = 4 ns and the analyzed number of frequencies was expanded to N ω = 2 k , k ∈ {3; 4; …; 9}, which gives f 0 = 62.5 MHz and f N ∈ {1; 2; 4; 8; 16; 32; 64} GHz.Fig. Accuracy of TPSF generation with TD -time domain and FD -frequency domain methods using low resolution "slab LO" and high resolution "slab HI" homogenous FEM models.Residuals are related to analytical solution of diffusion equation within the 80%-5% TPSF region of interest.N ω -number of frequencies of FD method, N t -number of time steps TD method, sd -source-detector separation.
Figure 5 shows that the frequency domain approximation error reduces much faster as compared to the time domain solution.Further, the low resolution model slab LO may be considered inaccurate for short source-detector separation r sd = mm as the field is rapidly changing near the sources and the model to be appropriately represented to ensure low numerical of Galerkin approximation at short source-detector separations.The frequency domain solution appears to stable for all r sd values, which is not the case for the time domain Table 3 shows calculation times per source of the spatial distributions of time-resolved photon fluence rate Φ(r,t) within the finite element models as in Fig. 2. The calculation time is linearly proportional to t and N ω in time and frequency domain respectively.shown in Table 3, the lowest calculation time within 5% error range in the frequency domain is ~50% than in the time domain.In other words, the same accuracy can be achieved ~2 faster when utilizing the proposed formulation based on frequency decomposition.difference in execution time between time and frequency domains for the same number of instances of solving the forward problem solver (where N t = N ω ) is a consequence of using complex numbers in the frequency domain, where arithmetic of complex numbers doubles the required execution time.4 shows errors of recovered optical parameters using slab LO model and both methods at N t = 128 time points and N ω = 32 frequencies respectively.Input data was generated with analytical solution of diffusion equation and source-detector separation r sd = 29 mm.Observation time T was calculated for each TPSF curve as 3 times the mean time of flight of photons (T = 3⟨t⟩).A shot noise was added to the data assuming that each TPSF curve is built with 10 7 photons.The nonlinear curve fitting algorithm (interiorreflective Newton method) used which always started μ a = 0.1 −1 , μ′ s = 2 mm −1 and set with error function and step tolerance of 10 −6 and the following constrains: μ a ∈ [0; 0.1] mm −1 , μ′ s ∈ [0.05; 2] mm −1 .Results show similar recovery errors for both methods.The fitting algorithm converged within the same number of iterations for both methods (5 to 10 depending on assumed optical properties).The heterogeneous finite element models in Fig. 2 consist of highly absorbing structures (veins, arteries) and low absorbing and scattering regions (trachea, cerebrospinal fluid).As a consequence the condition number of the forward problem (AΦ = q) described by N × N sparse matrix A [25] is high, where N is the number of finite element mesh nodes.However, in this work iterative solvers are deployed, using linear system preconditioning as follows: T T and , = = MAM y Mq Φ M y (6) where Φ N × 1 is the unknown photon fluence rate, q N × 1 is a light source vector and M N × N is factorized sparse preconditioner, which satisfies following conditions: MAM T ≈ I and M T M ≈ A −1 .Equation ( 6) is then solved using bi-conjugate gradient stabilized iterative solver [44].Table 5 compares execution time per source as a function of different number of mesh nodes.The condition numbers is the 1-norm estimate (MATLAB 'condest' function).As shown in Table 5 the execution time is also related to the condition number of the preconditioned system MAM T .The Factorized Sparse Approximate Inverse preconditioning method has been utilized in this work and modified to be used in massive parallel calculation environment.The parallel preconditioning is extremely fast (few milliseconds) when compared to the time needed to solve Eq. ( 6).However, it does not reduce the problem condition number linearly in relation to any of parameters presented in Table 5.

Discussion and conclusions
This work is motivated by a need of fast analysis of time-resolved data using realistic, heterogeneous tissue models.The time-resolved near infrared spectroscopic model built on the frequency domain approximation executed on a (GPU) accelerated platform opens a path to a reliable real-time curve fitting analysis of time-resolved data for heterogeneous tissue models.
The proposed method provides the most accurate result when the time T is linked to the expected region of interest, e.g. by using a priori knowledge of the source-detector distance rsd and/or scattering and absorption properties.It is proposed that time T can be set to T = 3⟨t⟩.The mean time of flight of photons ⟨t⟩ can be approximately assessed as ⟨t⟩ = DPF⋅r sd ⋅n/c 0 or more precisely ⟨t⟩ = r sd 2 /[2c 0 n −1 D 0.5 (r sd μ a 0.5 + D 0.5 )] [11] when optical properties can be assumed a priori.c 0 / mm⋅s −1 is the speed of light in vacuum, n / -is the medium refractive index, r sd / mm is investigated source-detector separation, μ a / mm −1 is the absorption coefficient, D / mm is the diffusion coefficient and the differential pathlength factor (DPF / -) as reported in [45,46] can be set to 6 for a human head.The DPF depends on for example age [47] and should be used with care.If a better representation of late photons fluence rate is needed, that is the 80%-5% region is changed to 80%-1% region or more, the number of frequencies N ω should be increased.The 80% limit cuts off the early photons where the diffusion approximation does not represent the light transport well.The late photons limit is often related to the noise level of measured curve and may vary between 20% -1% [40][41][42][43].
The proposed methodology highlights which frequency components of modulated or pulsed light sources carry usable information about tissue hemodynamic properties.The usable frequency band for tissue models as described in this work can be limited to [0.125; 8] GHz.It has been shown in [48] that for a multi-frequency diffuse optical tomography 13 frequencies between 20 MHz and 500 MHz can be used to improve tomographic reconstruction.However, considering for example a 30 mm source-detector separation and the neck model, frequencies lower than 125 MHz can carry redundant information as shown in Fig. 1.Low frequencies will not affect the shape of the TPSF significantly since cosines of

Fig. 1 .
Fig. 1.A schematic representation of contributions to a TPSF curve of angular frequencies close to ω = π/(2T) for maximum of a TPSF at t = 0 (a) and t = T/3 (b).T -maximum observable photons travel time, A -cosines amplitude.

Fig. 2 .
Fig. 2. Finite element models/meshes of homogenous slab (a), human head (b) and human neck (c).Green cubes show the 8 sources and red spheres show 9 detectors organized into a grid with 13 mm inter-optode distance.Models are not in scale.