A novel straightforward, accessible and efficient approach is presented for performing hyperspectral time-domain diffuse optical spectroscopy to determine the optical properties of samples accurately using geometry specific models. To allow bulk parameter recovery from measured spectra, a set of libraries based on a numerical model of the domain being investigated is developed as opposed to the conventional approach of using an analytical semi-infinite slab approximation, which is known and shown to introduce boundary effects. Results demonstrate that the method improves the accuracy of derived spectrally varying optical properties over the use of the semi-infinite approximation.
© 2016 Optical Society of America
Diffuse Optical Spectroscopy (DOS) is becoming a routine tool for the optical characterization of a variety of diffusive materials, allowing the analysis of biological tissues , tissue phantoms , wood [3–5], food and powders [6,7]. In a biological application DOS systems provide useful information regarding the pathophysiological status of tissue for both diagnostic and therapeutic applications. Recently, emphasis has been on the development and application of instruments that can non-invasively provide the spectral absorption coefficient, μa, and the reduced scattering coefficient, μs’, for near-infrared (NIR) wavelengths (~600–1000nm). Multi-wavelength DOS systems are preferential to single wavelength systems as they yield more robust and accurate measurements, particularly in the presence of noise . Spectroscopic approaches also allow the quantification of chromophores based on their known absorption spectra. The four major contributors to absorption in biological tissues are deoxyhaemoglobin, oxyhaemoglobin, lipids and water, and spectroscopic methods can also cater to other tissue chromophores [9,10].
Time-resolved (TR-) DOS is one technique for spectroscopically measuring bulk μa and μs’, typically across a broad range of NIR wavelengths. It involves the use of ultrafast laser pulses and single-photon time-resolved detection to measure the temporal point spread function (TPSF) of transmitted or diffusely reflected light. Optical properties are then calculated by minimising the difference between measured TPSFs and those obtained by a model of light transport through the sample. Conventionally TPSFs are modelled analytically by assuming a semi-infinite slab geometry. However, this is often a poor approximation and errors are encountered due to unexpected losses at boundaries (so-called boundary effects) . This limits accuracy and makes results dependent upon the particular experimental configuration as different boundary effects are encountered. In some cases, this becomes especially severe such as when samples are small and/or highly irregular (e.g. small animal imaging phantoms).
Figure 1 shows an illustration of boundary effects for a 2D circular sample modelled as a semi-infinite slab. The figure shows the simulated absorption sensitivity for the single source-detector configuration, calculated using NIRFAST  assuming homogeneous μa = 0.01mm−1 and μs’ = 1mm−1. The colours indicate the (log-scale) sensitivity to spatially resolved changes in μa and plotted contour lines effectively indicate the average photon path through the slab model (with decreasing number of photons per contour), the dotted line shows a circular geometry for which light is lost (through the boundary) relative to the semi-infinite model.
Whilst the effects in the case shown in Fig. 1 may appear to be small, the light measurement made at the detector integrates over all boundary effects so the impact is still significant. Moreover, when operating in the time domain, long-lived photons that travel further in the medium are more affected by boundary effects than early photons, causing a distortion on the TPSF. To provide an indication of the size of the effects, Fig. 2 shows a simple example in the form of the relative magnitude of the simulated (continuous wave) measurement for the geometry of Fig. 1 as compared to that of a semi-infinite slab with identical optode configuration. This therefore quantifies the boundary loss in amplitude as a function of bulk optical absorption and reduced scattering coefficients.
From Fig. 2, it can be seen that the impact of the boundary effects is higher with decreasing absorption and that for low (but still realistic) levels of μa (for biological tissues in NIR and for tissue phantoms) it is greater than 20%. It intuitively follows that the boundary effects are less relevant in cases of higher absorption because less light is able to reach the boundary in a more highly attenuating case (in this geometry the photon distribution is more confined since strong absorption will kill the longest photon paths). The relationship with increasing μs' is more complicated and non-linear owing to the more ballistic behaviour at low scattering (in the present geometry this reduces the relevance of the boundaries) coupled with the attenuating effect of high scattering. For the broad range of scattering values shown here, the impact of the boundary effects is relatively constant at 16-22%. It therefore appears that the boundary effect impact is more sensitive to absorption than scattering.
The most obvious solution to overcome this issue and grant accurate estimation of bulk properties in finite-size samples would be to use an accurate shape-specific forward solver (either using a finite element model (FEM) or Monte-Carlo (MC) based approach) and a model-based inversion procedure . However, this approach can be resource-intensive, both computationally and in terms of time, as it requires the construction of complex inverse problem software and involves the generation of data for many combinations of μa and μs for each measurement point. In the case of spectral measurements or experiments involving many samples this can be inefficient.
The goal of this paper is to present a novel straightforward, accessible and efficient method for the recovery of bulk parameters from TR-DOS measurements of arbitrary geometries. To speed-up and facilitate the process, we adopt a method initially proposed for fitting with MC simulations that was previously demonstrated solely for a simple slab geometry [14–16]. The idea is to generate (only once) for each specific geometry a finite set of purely scattering simulations that can be easily exploited in a fitting algorithm to quickly reproduce the effect of any scattering coefficient by linear interpolation of the library curves and any absorption coefficient by multiplication for the Lambert-Beer factor. Models are built as finite element meshes within which the diffusion equation (DE) is used to model light transport. The advantages of the presented technique are that the model library can be built on-the-fly very quickly as compared to using Monte Carlo, and with significantly better geometric accuracy (accounting for boundary effects) compared to simple analytical models. As a consequence, it is found that fitting results using this method are improved for small and irregular samples.
2. Materials and methods
2.1. Tissue phantoms
Three solid polyurethane tissue phantoms were studied in simulation and practical experiments (Fig. 3). These were a block phantom (INO, Canada), a cylinder phantom (INO) and a mouse phantom  (XPM2, Perkin Elmer), each having nominally spatially homogeneous, spectrally varying optical properties similar in magnitude to those of biological tissues . The mouse phantom in particular has a wavelength-dependent absorption spectrum similar to that of blood. Absolute spectrally varying absorption and scatter values are known for the block and circular phantoms, and a (low-resolution) graphical representation of these is available for the mouse phantom, from which absorption and scatter values can be manually determined. For small samples such as these, a transmittance measurement geometry is often preferable to a reflectance scheme due to the natural suppression of stray light between the injection and detection fibres and easier placement of the probes, particularly in the case of irregularly shaped specimens. For these reasons, we investigated the transmittance-mode geometries shown in Fig. 3 both for simulations and experiments. For extra insight, in experiment we carried out further measurement in reflectance mode (for the mouse phantom only as shown in Fig. 3(c)).
2.2. Geometry-specific finite element modelling
Physical phantoms were modelled using geometrically accurate FEM meshes. In the case of the block (43631 nodes corresponding to 232913 linear tetrahedral elements) and cylinder (7530 nodes and 30746 elements) phantoms, meshes were easily generated due to the known simple geometry. For the mouse phantom, a surface mesh provided by the manufacturer was used from which a volumetric FEM (35058 nodes and 176514 elements) was generated. Source and detector positions were modelled using measurements of physical fibre positions made by hand. These were at 26, 25 and 22 mm separation distances for the block, cylinder and XPM2 phantom respectively.
Open-source FEM software (NIRFAST ) was used for geometry-specific light modelling. Synthetic forward data was created using a single set of spectral optical properties for simulation studies (see Fig. 5 for nominal values) and a TPSF library was constructed for each phantom for use in the fitting routines for simulation and experimental data. TPSFs were calculated using time-resolved methods based on the diffusion approximation . Each library comprised 21 models for μs' ranging from 0.1 to 5mm−1 (in steps of 0.245mm−1) and μa = 0mm−1 for a total of 6ns with time-bins of 10 ps. The total computation time for these sets of TPSFs was ~600 s, on a MacBook Pro with a 2.9 GHz Intel Core i7, 8 GB 1600 MHz DDR3 RAM running MATLAB 2015.
2.3. Time-resolved diffuse optical spectroscopy (TR-DOS)
TR-DOS experiments were performed using a previously described system [20,21] (Fig. 4). Briefly, the system comprises a broadband, pulsed laser source coupled to a rotating prism acting as a tunable bandpass filter. Spectrally filtered light is injected into the sample by an optical fiber in contact with its surface. Detection is via another contact-mode fiber, which is coupled directly into the detector. Time-resolved detection is performed by a Hybrid PMT (HPM-100-50, Becker and Hickl, Germany) providing an efficiency of ~15% for wavelengths in the range 400-820nm, and a temporal resolution <150 ps.
To take into account the temporal width of the laser source and the dispersion caused by the optical fibers and detectors, an Instrumental Response Function (IRF) was acquired by directly coupling the detection and collection fibers.
2.4. Optical property fitting procedure
The optical property fitting algorithm is based on a Levenberg-Marquardt minimization  of experimental (or simulated) and model TPSFs convolved by the measured IRF. Forward model library data (for a discrete set of scattering coefficients with μa = 0) is used to optimize μs' (missing scattering values are considered via a linear interpolation) and absorption is added analytically using multiplication by the factor , where v is the light speed inside the sample and t is the photon arrival time. This approach has previously been demonstrated using TPSFs generated from a library of Monte Carlo solutions . In this study the fitting is performed firstly based on semi-infinite slab model data and secondly using geometrically accurate numerical diffusion model data. The part of the TPSF used in all fittings (practical and simulation) extended from 80% of the peak on the leading edge, down to 1% on the falling edge.
Three simulation experiments were conducted featuring each of the phantoms shown in Fig. 3, in transmittance-mode only with source-detector configurations as outlined above. For simplicity, identical optical properties were assigned to each phantom. These were calculated to represent a biological tissue containing 0.01 mM oxyhaemoglobin, 0.01 mM deoxyhaemoglobin and 40% water, with spectral power-law scattering of the form μs’ = Aλ-b for A = B = 1, with λ expressed in micrometres and μs’ in mm−1. Forward simulated data was generated using NIRFAST and TR-DOS fitting was applied using a model-specific NIRFAST-generated TPSF library as well as the conventional semi-infinite slab library.
Figure 5 shows fitted optical properties obtained for each phantom using each of the models, along with true target values. There is good agreement between the fitted data and the true values in all cases where the geometrically accurate model was used. Where the semi-infinite approximation was applied to the block phantom there is no error on μs’ but an overestimation of μa. The finite lateral dimensions of the phantom (16.5 mm from the fiber axis) can be considered to account for photon losses at the boundaries. For the cylinder and the mouse phantoms there is a larger level of error in μa, with recovered values being approximately double the true values. In μs’ the mouse phantom data shows slightly higher error than the cylinder with the minimum error being approximately 8% as opposed to 4%. There is also a strong distortion of the reduced scattering spectrum coinciding with the spectral location of the absorption peak which may be due to the complex interplay between absorption, scattering and boundary effects outlined in section 1 and Fig. 2. Of the observed errors, those in μs' are relatively small compared to those in μa. This is consistent with the expectation that μa has a higher sensitivity to boundary effects (Fig. 2).
Overall, it appears that as the complexity in the phantom geometry is increased, the semi-infinite model becomes less applicable and the results are more severely distorted owing to boundary effects. The geometrically accurate model approach appears to effectively remove the influence of these and return fitting results that are geometry-independent and correct.
Figure 6 shows recovered optical properties obtained using the semi-infinite slab and geometrically accurate models for practical experiments performed using the block and cylinder phantoms. For the block phantom there are no major differences between the results obtained using the geometrically accurate model and those for the semi-infinite approximation; across all wavelengths, μa and μs’ results are effectively identical. These also match reasonably well with the nominal values for the physical phantom, though the match is less accurate in μa at shorter wavelengths where the maximum difference is approximately 20%. Some erratic values are observed in the nominal μs’, which are assumed to be noise effects owing to the fact that scattering spectra are typically smooth and non-increasing and the dynamic range is small.
For the cylinder phantom, the μa values recovered with the semi-infinite model are approximately 40% higher than the manufacturer’s specification. With the geometrically accurate model the values are closer to the manufacturer data, to within approximately 10% across all wavelengths, demonstrating a substantial improvement. The fact that μa is reduced is consistent with expectations related to boundary effects (as seen in simulation experiments) as these artificially raise the calculated μa when using the semi-infinite model. The recovered scattering values for the cylinder are approximately 3% higher than the manufacturer data and are increased by a further 3% when using the geometrically accurate model, although this is a relatively insignificant change.
Figure 7 shows results for the mouse phantom for the two experimental configurations shown in Fig. 3(c); a transmittance measurement through the centre of the phantom and a reflectance measurement with source-detector separation distance of 10mm. Owing to the evident discrepancy between all fitting results and the manufacturer-provided nominal data, three additional “validation” data sets are included. These are repeat measurements made using the standard (semi-infinite) fitting method on reflectance data with very short source-detector separation distance (5mm). Owing to the short path length in these cases any boundary effects are expected to be minimal and as such these values are expected to be a highly accurate depiction of the ground truth.
For the transmittance-mode data, μa values obtained using the semi-infinite model are slightly below the manufacturer data and significantly reduced when using the geometrically accurate model. This is consistent with the simulation and other data in that the recovered magnitude of the values are reduced when applying the geometrically accurate model (consistent with boundary effect elimination) though it is surprising that the nominal values are higher than both sets of the fitted data. The validation measurements are in good agreement with one another and those of the geometrically-accurate model, suggesting that the nominal values may be at fault. Fitted scattering values vary little between the two fitting approaches but are all approximately 50% higher than the nominal values, however the validation data is higher still by approximately the same margin. Again this indicates an error in the manufacturer data or its interpretation (given that it was manually extracted from a small chart).
For the reflectance-mode experiment, the geometrically-accurate model absorption results (Fig. 7(b)) are near-identical to those obtained in transmittance-mode. They are therefore also in good agreement with the validation data. This is encouraging as it shows independence of the absorption measurement with respect to the experimental configuration, due to boundary effect elimination. Once again the semi-infinite model values are higher, though not as high as in the transmittance-mode case. This is as expected because of the shorter source-detector distance and reflectance geometry making the semi-infinite approximation more valid in this case and as such less (though still somewhat) susceptible to boundary effects. In scattering, both the semi-infinite and geometrically-accurate approaches return similar values, as before, although these are in better agreement with the validation data than the transmittance-mode counterparts. One possible reason for this is that the mouse phantom actually contains an internal (deactivated) light source  creating an unknown optical anomaly, it is possible that this may have interfered with the scattering estimate, which could explain why both semi-infinite and geometry matched models returned near-identical results but both reduced in transmittance-mode. Another possibility is that these are boundary effects that are harder to correct; influenced by μs’ being less sensitive to boundary effects (Fig. 2).
4. Discussions and conclusions
A novel practical approach for performing time-resolved diffuse optical spectroscopy (TR-DOS) using geometrically accurate light models has been investigated with the aim of eliminating the detrimental influence of boundary effects on fitted bulk optical properties (μa and μs’) at little added cost. The proposed method incorporates subject-specific temporal point spread functions (TPSFs) using geometrically accurate finite element models as opposed to the simple but geometrically inaccurate semi-infinite slab models typically employed. This was achieved using open source light modelling software [12,23] which uses the diffusion approximation to make fast calculations in a practical time-frame: In the work presented here, for the mouse phantom, the total time for creating the model specific TPSF libraries was approximately 10 minutes using a laptop with a 2.9 GHz Intel Core i7 and 8 GB of RAM.
Three biological tissue mimicking phantoms with realistic optical properties were studied in simulation and experiments. Simulations (Fig. 5) showed that as the sample geometry became increasingly complex, the recovered bulk optical parameters obtained with a semi-infinite slab model were increasingly influenced by boundary effects, leading to greater errors in recovered μa and μs’. The presented approach appeared to remove these effects, providing a more accurate match to the true values.
Using experimental data for a range of phantoms, similar results were obtained for the block and cylindrical phantoms. The more complex mouse shaped phantom was studied in two experimental configurations; transmittance and reflectance-mode. Across the two setups the semi-infinite model results differed due to (different) boundary effects whereas the recovered μa values obtained with the new method were consistent. Although the results did not match expected values, validation data agreed with the fitting, indicating erroneous nominal data. Across all three experiments, μs' showed very little change between the semi-infinite and geometrically-accurate models due to the fact that μs’ is less sensitive to boundary effects (Fig. 2).
In general, presented approach has been demonstrated to improve the accuracy of TR-DOS when working with homogeneous samples of complex geometry. The only requirement is for a geometrically accurate numerical depiction of the geometry. Due to the widespread availability of 3D scanning systems, coupled with the availability of open source meshing software (e.g. NIRVIEW ), this is now readily achievable and the approach can find immediate application. One possible application is the characterization of tissue-like phantoms as demonstrated in this work, this approach may be more accurate and time-efficient than making measurements on large slabs of excess material as is otherwise required .
Considering the application of this method to heterogeneous biological samples, it is useful to mention that the use of a homogenous model does not necessarily lead to an “average” estimate of the heterogeneous sample properties. The type of heterogeneity (layered, localized), measurement configuration (reflectance, transmittance) and optical properties will affect the outcome in combination. This is an existing issue for current TR-DOS approaches (semi-infinite or slab specimens) and no more so here. However, the proposed approach could potentially be extended to deal with segmented images of heterogeneous media. In this case, the identification of a small number of regions would permit the construction of a library with a limited number of scattering combinations (for each region) to be used in a fast reconstruction scheme where absorption properties are accounted applying the Lambert-Beer factor linked to the mean time spent in each region. This approach deserves specific studies in the future.
This work has been funded by Laserlab-Europe (EU-FP7 284464) and EPSRC (EP/F50053X/1).
References and links
2. 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] [PubMed]
3. R. Kitamura, T. Inagaki, and S. Tsuchikawa, “Determination of true optical absorption and scattering coefficient of wooden cell wall substance by time-of-flight near infrared spectroscopy,” Opt. Express 24, 569–574 (2016).
4. A. Farina, I. Bargigia, E.-R. Janeček, Z. Walsh, C. D’Andrea, A. Nevin, M. Ramage, O. A. Scherman, and A. Pifferi, “Nondestructive optical detection of monomer uptake in wood polymer composites,” Opt. Lett. 39(2), 228–231 (2014). [CrossRef] [PubMed]
5. S. Tsuchikawa and M. Schwanninger, “A review of recent near infrared research for wood and paper (Part 2),” Appl. Spectrosc. Rev. 48(7), 560–587 (2013). [CrossRef]
6. R. Cubeddu, C. D’Andrea, A. Pifferi, P. Taroni, A. Torricelli, G. Valentini, C. Dover, D. Johnson, M. Ruiz-Altisent, and C. Valero, “Nondestructive quantification of chemical and physical properties of fruits by time-resolved reflectance spectroscopy in the wavelength range 650-1000 nm,” Appl. Opt. 40(4), 538–543 (2001). [CrossRef] [PubMed]
7. C. Abrahamsson, J. Johansson, S. Andersson-Engels, S. Svanberg, and S. Folestad, “Time-resolved NIR spectroscopy for quantitative analysis of intact pharmaceutical tablets,” Anal. Chem. 77(4), 1055–1059 (2005). [CrossRef] [PubMed]
8. C. D’Andrea, L. Spinelli, A. Bassi, A. Giusto, D. Contini, J. Swartling, A. Torricelli, and R. Cubeddu, “Time-resolved spectrally constrained method for the quantification of chromophore concentrations and scattering parameters in diffusing media,” Opt. Express 14(5), 1888–1898 (2006). [CrossRef] [PubMed]
9. 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] [PubMed]
10. P. Taroni, D. Comelli, A. Pifferi, A. Torricelli, and R. Cubeddu, “Absorption of collagen: effects on the estimate of breast composition and related diagnostic implications,” J. Biomed. Opt. 12(1), 014021 (2007). [CrossRef] [PubMed]
12. 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. Meth. Engng. 25, 711–732 (2009).
13. G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt. 45(5), 1062–1071 (2006). [CrossRef] [PubMed]
16. 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). [CrossRef] [PubMed]
17. C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging,” J. Biomed. Opt. 12(2), 024007 (2007). [CrossRef] [PubMed]
18. W. F. Cheong, S. A. S. Prahl, and A. J. A. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. 26(12), 2166–2185 (1990). [CrossRef]
19. 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). [CrossRef] [PubMed]
20. A. Bassi, A. Farina, C. D’Andrea, A. Pifferi, G. Valentini, and R. Cubeddu, “Portable, large-bandwidth time-resolved system for diffuse optical spectroscopy,” Opt. Express 15(22), 14482–14487 (2007). [CrossRef] [PubMed]
21. I. Bargigia, A. Tosi, A. Bahgat Shehata, A. Della Frera, A. Farina, A. Bassi, P. Taroni, A. Dalla Mora, F. Zappa, R. Cubeddu, and A. Pifferi, “Time-resolved diffuse optical spectroscopy up to 1700 nm by means of a time-gated InGaAs/InP single-photon avalanche diode,” Appl. Spectrosc. 66(8), 944–950 (2012). [CrossRef] [PubMed]
22. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 1988).
23. 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). [CrossRef] [PubMed]