Open Access
Add to BibSonomy
Abstract
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
1. Introduction
Time-resolved spectroscopy (TRS) devices [1–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 two-homogeneous 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 time-resolved photon fluence rate within MRI based head model (1 mm3 voxel size) in about 1.7 min for 107 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.
2. 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.
2.1 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 r is a point in R3 Euclidean space, D(r) = [3(μa(r) + μ′s(r))]−1 is the diffusion coefficient, μa(r) is absorption coefficient, μ′s(r) is the reduced scattering coefficient, t is time, Φ(r,t) is photon fluence rate (TPSF), cm(r) is speed of light in the medium which is a function of the medium’s refractive index n(r), rs is source location and S(rs,t) is a source term e.g. a pulse of light (Dirac delta δ(rs,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 Nt points, where the source term S(rs,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 F of the source term δ(rs,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ω.2.2 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).
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.
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 rsd > 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−2s−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.
3. 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 AN × N is a sparse matrix representing finite element model defined by N vertex nodes in R3 Euclidean space, qN × 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.
3.1 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.
Table 1. Quality measures of finite element models.
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.
Table 2. Optical properties of the neck and head models at 780 nm [24, 33, 34].
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.
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.
3.2 Accuracy and speed
Numerical solvers for the time-resolved forward problem as defined in Eq. (2) use a solution of the previous time points tn-s (s ≥ 1) to solve the problem at point tn. Thus, the solution accuracy increases for shorter integration steps Δt = tn - tn-1 = T/Nt (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 Nt = 2k, 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 Nt = 210. 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–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 Nt of NIRFAST time domain method as defined in Eq. (2). Data shown for the neck model as in Fig. 2c and two source-detector pairs. rsd – 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 Nt = 1024. Residuals are related to TPSFs calculated for Nt = 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ω = 2k, 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 Nt = 1024. The lowest usable frequency needed to reproduce the TPSF is set to f0 = ω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: fN = ω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. 2c and two source-detector pairs. rsd – 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 Nt = 1024 Residuals are related to TPSFs calculated in time domain for Nt = 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 rsd 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ω = 2k, k ∈ {3; 4; …; 9}, which gives f0 = 62.5 MHz and fN ∈ {1; 2; 4; 8; 16; 32; 64} GHz.
Fig. 5 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, Nt – number of time steps of TD method, rsd – 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 as inaccurate for the short source-detector separation rsd = 13 mm as the field is rapidly changing near the sources and the model needs to be appropriately represented to ensure low numerical errors of Galerkin approximation at short source-detector separations. The frequency domain solution appears to be stable for all rsd values, which is not the case for the time domain solution.
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 shown in Fig. 2. The calculation time is linearly proportional to Nt and Nω in time and frequency domain respectively. As shown in Table 3, the lowest calculation time within the ± 5% error range in the frequency domain is ~50% lower than in the time domain. In other words, the same accuracy can be achieved ~2 times faster when utilizing the proposed formulation based on frequency decomposition. The difference in execution time between time and frequency domains for the same number of instances of solving the forward problem solver (where Nt = Nω) is a consequence of using complex numbers in the frequency domain, where arithmetic of complex numbers doubles the required execution time.
Table 3. Calculation times per source of 3D distributions of time-resolved photon fluence rate Φ(r,t) within the finite element models. Marked cells show the lowest calculation time within ± 5% relative error. TD – time domain, FD – frequency domain, Nt – number of time point in TD and Nω – number of frequencies in FD.
Table 4 shows errors of recovered optical parameters using slab LO model and both methods at Nt = 128 time points and Nω = 32 frequencies respectively. Input data was generated with analytical solution of diffusion equation and source-detector separation rsd = 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 107 photons. The nonlinear curve fitting algorithm (interior-reflective Newton method) was used which always started at μa = 0.1 mm−1, μ′s = 2 mm−1 and was set with error function tolerance 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).
Table 4. Error of optical parameters recovery using slab models and TD – time domain, FD – frequency domain methods at Nt = 128 time point and Nω = 32 frequencies respectively.
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:
where ΦN × 1 is the unknown photon fluence rate, qN × 1 is a light source vector and MN × N is factorized sparse preconditioner, which satisfies following conditions: MAMT ≈ I and MTM ≈ 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 MAMT. 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.Table 5. Execution time per source per mesh node for the finite element models. TD – time domain, FD – frequency domain.
4. 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⋅rsd⋅n/c0 or more precisely ⟨t⟩ = rsd2/[2c0n−1D0.5(rsdμa0.5 + D0.5)] [11] when optical properties can be assumed a priori. c0 / mm⋅s−1 is the speed of light in vacuum, n / - is the medium refractive index, rsd / 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–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 low frequencies will be almost flat curves within the observation window [0; T]. The upper frequency limit ωN = 2Nωω0 is defined by the minimum number of frequencies Nω = 32 required to reproduce the TPSF. The frequency band utilized will change with source-detector separation and tissue scattering. Other factors, which influence the TPSF (including tissue absorption), will have lower impact on the frequency band.
Theoretically, it is possible to recover a distribution of time of flight of photons using a frequency domain instrument capable of measurements on 32 frequencies within [0.125; 8] GHz frequency band. It has been discussed in [49] that the frequency domain measurement cannot easily differentiate between long pathlength photons, which is understandable since the time of flight of a photon through a tissue should be lower than the period of a modulated source to avoid phase overlap (phase shift ≥ 2π). However, this restriction does not influence the proposed superposition of frequency components which is independent of the phase shift overlap.
All calculations were carried out with graphical processing unit (GPU) accelerated NIRFAST version. Speed of GPU parallel solvers is highly dependent on hardware. The execution time shown in this paper can be reduced using more efficient or multiple GPUs.
Funding
European Union‘s Horizon 2020 research and innovation programme (88303, LUCA project); Fundació CELLEX Barcelona; the “Severo Ochoa” Programme for Centres of Excellence in R&D (SEV-2015-0522); the Obra social “la Caixa” Foundation (LlumMedBcn) and AGAUR-Generalitat (2014SGR-1555).
Disclosures
The authors declare that there are no conflicts of interest related to this article. LUCA project involves industrial collaboration and, as such, potential conflicts of interest are being monitoring by relevant institutional bodies. None has been identified to date.
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). [PubMed]
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). [PubMed]
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). [PubMed]
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). [PubMed]
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). [PubMed]
6. D. Grosenick, H. Rinneberg, R. Cubeddu, and P. Taroni, “Review of optical breast imaging and spectroscopy,” J. Biomed. Opt. 21(9), 091311 (2016). [PubMed]
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). [PubMed]
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). [PubMed]
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). [PubMed]
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). [PubMed]
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). [PubMed]
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). [PubMed]
14. A. Puszka, L. Hervé, A. Planat-Chrétien, A. Koenig, J. Derouard, and J.-M. Dinten, “Time-domain reflectance diffuse optical tomography with Mellin-Laplace transform for experimental detection and depth localization of a single absorbing inclusion,” Biomed. Opt. Express 4(4), 569–583 (2013). [PubMed]
15. 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). [PubMed]
16. J. A. Guggenheim, I. Bargigia, A. Farina, A. Pifferi, and H. Dehghani, “Time resolved diffuse optical spectroscopy with geometrically accurate models for bulk parameter recovery,” Biomed. Opt. Express 7(9), 3784–3794 (2016). [PubMed]
17. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 50902 (2013). [PubMed]
18. L. Wang, S. L. Jacques, and L. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [PubMed]
19. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [PubMed]
20. H. Dehghani, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Optical tomography in the presence of void regions,” J. Opt. Soc. Am. A 17(9), 1659–1670 (2000). [PubMed]
21. Q. Fang and D. A. Boas, “Monte Carlo Simulation of Photon Migration in 3D Turbid Media Accelerated by Graphics Processing Units,” Opt. Express 17(22), 20178–20190 (2009). [PubMed]
22. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20(2 Pt 1), 299–309 (1993). [PubMed]
23. 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). [PubMed]
24. M. Jermyn, H. Ghadyani, M. A. Mastanduno, W. Turner, S. C. Davis, H. Dehghani, and B. W. Pogue, “Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography,” J. Biomed. Opt. 18(8), 086007 (2013). [PubMed]
25. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2008). [PubMed]
26. Q. Zhu, H. Dehghani, K. M. Tichauer, R. W. Holt, K. Vishwanath, F. Leblond, and B. W. Pogue, “A three-dimensional finite element model and image reconstruction algorithm for time-domain fluorescence imaging in highly scattering media,” Phys. Med. Biol. 56(23), 7419–7434 (2011). [PubMed]
27. A. Kienle, T. Glanzmann, G. Wagnières, and H. Bergh, “Investigation of two-layered turbid media with time-resolved reflectance,” Appl. Opt. 37(28), 6852–6862 (1998). [PubMed]
28. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnivres, and H. van den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37(4), 779–791 (1998). [PubMed]
29. M. Kacprzak, A. Liebert, P. Sawosz, N. Żolek, and R. Maniewski, “Time-resolved optical imager for assessment of cerebral oxygenation,” J. Biomed. Opt. 12, 034019 (2007).
30. A. Liebert, P. Sawosz, D. Milej, M. Kacprzak, W. Weigl, M. Botwicz, J. Mączewska, K. Fronczewska, E. Mayzner-Zawadzka, L. Królicki, and R. Maniewski, “Assessment of inflow and washout of indocyanine green in the adult human brain by monitoring of diffuse reflectance at large source-detector separation,” J. Biomed. Opt. 16, 046011 (2011).
31. J. Ashburner and K. J. Friston, “Image segmentation,” in Human Brain Function, 2 ed., R. S. J. Frackowiak, K. J. Friston, C. Frith, R. Dolan, K. J. Friston, C. J. Price, S. Zeki, J. Ashburner, and W. D. Penny, eds. (Academic Press, 2003).
32. K. J. Friston, “Introduction: experimental design and statistical parametric mapping,” in Human Brain Function, 2 ed., R. S. J. Frackowiak, K. J. Friston, C. Frith, R. Dolan, K. J. Friston, C. J. Price, S. Zeki, J. Ashburner, and W. D. Penny, eds. (Academic Press, 2003).
33. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [PubMed]
34. H. Fujii, Y. Yamada, K. Kobayashi, M. Watanabe, and Y. Hoshi, “Modeling of light propagation in the human neck for diagnoses of thyroid cancers by diffuse optical tomography,” Int. J. Numer. Methods Biomed. Eng. 33(5), 2826 (2016). [PubMed]
35. “CGAL, Computational Geometry Algorithms Library”, retrieved http://www.cgal.org.
36. A. Custo, W. M. Wells 3rd, A. H. Barnett, E. M. C. Hillman, and D. A. Boas, “Effective scattering coefficient of the cerebral spinal fluid in adult head models for diffuse optical imaging,” Appl. Opt. 45(19), 4747–4755 (2006). [PubMed]
37. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [PubMed]
38. G. Strangman, M. A. Franceschini, and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” Neuroimage 18(4), 865–879 (2003). [PubMed]
39. A. H. Hielscher, S. L. Jacques, L. Wang, and F. K. Tittel, “The influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissues,” Phys. Med. Biol. 40(11), 1957–1975 (1995). [PubMed]
40. M. Diop and K. St Lawrence, “Improving the depth sensitivity of time-resolved measurements by extracting the distribution of times-of-flight,” Biomed. Opt. Express 4(3), 447–459 (2013). [PubMed]
41. 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), 16010 (2014). [PubMed]
42. F. Martelli, S. Del Bianco, and G. Zaccanti, “Retrieval procedure for time-resolved near-infrared tissue spectroscopy based on the optimal estimation method,” Phys. Med. Biol. 57(10), 2915–2929 (2012). [PubMed]
43. 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). [PubMed]
44. H. A. d. Vorst, “Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems,” SIAM J. Sci. Statist. Comput. 13, 631–644 (1992).
45. A. Duncan, J. H. Meek, M. Clemence, C. E. Elwell, L. Tyszczuk, M. Cope, and D. T. Delpy, “Optical pathlength measurements on adult head, calf and forearm and the head of the newborn infant using phase resolved optical spectroscopy,” Phys. Med. Biol. 40(2), 295–304 (1995). [PubMed]
46. P. van der Zee, M. Cope, S. R. Arridge, M. Essenpreis, L. A. Potter, A. D. Edwards, J. S. Wyatt, D. C. McCormick, S. C. Roth, E. O. Reynolds, et al.., “Experimentally measured optical pathlengths for the adult head, calf and forearm and the head of the newborn infant as a function of inter optode spacing,” Adv. Exp. Med. Biol. 316, 143–153 (1992). [PubMed]
47. F. Scholkmann and M. Wolf, General Equation for the Differential Pathlength Factor of the Frontal Human Head Depending on Wavelength and Age (SPIE, 2013), 6.
48. X. Intes and B. Chance, “Multi-frequency diffuse optical tomography,” J. Mod. Opt. 52, 2139–2159 (2005).
49. T. Binzoni, A. Sassaroli, A. Torricelli, L. Spinelli, A. Farina, T. Durduran, S. Cavalieri, A. Pifferi, and F. Martelli, “Depth sensitivity of frequency domain optical measurements in diffusive media,” Biomed. Opt. Express 8(6), 2990–3004 (2017). [PubMed]
