1. Introduction

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 [1], tissue phantoms [2], 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 [8]. 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) [11]. 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 [12] assuming homogeneous μ_a = 0.01mm⁻¹ and μ_s’ = 1mm⁻¹. 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.

 

Fig. 1 Illustration of boundary effects for a circular 2D model, when modelled as a semi-infinite slab.

Download Full Size | PDF

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.

 

Fig. 2 Quantitative example of the impact of boundary effects on continuous wave measurements between a circular and semi-infinite model (a) as a function of μ_a and (b) as a function of μ_s', for a baseline value of μ_s' = 2mm⁻¹ and μ_a = 0.01mm⁻¹.

Download Full Size | PDF

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 [13]. 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 [17] (XPM2, Perkin Elmer), each having nominally spatially homogeneous, spectrally varying optical properties similar in magnitude to those of biological tissues [18]. 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)).

 

Fig. 3 Diagrams of tissue phantoms. (a) block phantom, (b) cylindrical phantom, (c) mouse-shaped phantom. The single studied transmittance-mode source-detector configuration is shown in each case along with the reflectance-mode configuration investigated for the mouse phantom only (marked “R”).

Download Full Size | PDF

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 [12]) 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 [19]. Each library comprised 21 models for μ_s' ranging from 0.1 to 5mm⁻¹ (in steps of 0.245mm⁻¹) and μ_a = 0mm⁻¹ 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.

 

Fig. 4 TR-DOS system schematic. Supercontinuum light is spectrally dispersed by a Pellin-Broca prism and then focused with an achromatic lens (f = 150 mm) on the tip of a graded-index optical fiber (50 μm core diameter). The fiber core acts as a spatial filter selecting a bandwidth of less than 10 nm. Transmitted light from the sample is collected by a step-index fiber (1 mm core diameter) and coupled to a Hybrid PMT connected to a Time-Correlated Single Photon Counting (TCSPC) board.

Download Full Size | PDF

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 [22] 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 $e^{-{\mu }_{a}vt}$, 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 [16]. 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.

3. Results

3.1. Simulation

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⁻¹. 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).

 

Fig. 5 Simulation results. (a-c), fitted absorption coefficient. (d-f), fitted reduced scattering coefficient. Results are shown for each phantom studied, for the geometry-specific and the semi-infinite slab based fitting approaches, along with the true values used to generate the forward data.

Download Full Size | PDF

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.

3.2. Experimental

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.

 

Fig. 6 Fitted absorption and scattering for the physical block and cylinder phantoms. Plots show nominal values and those obtained using the standard semi-infinite approximation and presented geometry-specific method as indicated in the legends.

Download Full Size | PDF

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.

 

Fig. 7 Fitted absorption and scattering for the physical mouse-shaped phantom in transmittance and reflectance configurations (R = reflectance mode, source-detector distance of 10mm). Plots show nominal values and those obtained using the standard semi-infinite approximation and presented geometry-specific method, with three additional validation data sets, as indicated in the legends. The validation data sets are three repeats of standard semi-infinite approximation based reflectance-mode measurements with a very short (5mm) source-detector separation (i.e. with minimal boundary effects).

Download Full Size | PDF

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 [17] 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 [23]), 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 [17].

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.