1. Introduction

In recent years, the possibility to characterize in-vivo and non-invasively the optical properties (absorption and scattering) of biological samples has attracted a great interest in the field of medical imaging. Pathologies like breast cancer [1,2] or osteoarticular diseases [3] are related to localized changes of optical properties due to increase of vascularization or collagen content.

Another important field is brain imaging. In this case, a map of oxy- and deoxy-hemoglobin concentration is fundamental for the diagnosis of injuries like ischemia [4], hemorrhage [5] and for functional imaging during a variety of tasks [6]. Recently, we have also proposed the monitoring of autologous tissues (“flap”) in reconstruction surgery using diffuse optical techniques. The vascularization of these tissues is fundamental and a fast postoperative control is important for the success of the operation [7].

A common approach to obtain a full 3D map of the optical properties in biological media at depths of a few cm is based on the combination of a set of light sources and detectors placed on the surface of the medium. By measuring the light collected at the detector d due to the light injected at the source s, it is possible by the help of sensitivity matrices to recover the distribution of optical properties inside the sample. This technique is called Diffuse Optical Tomography (DOT) [8].

Different modalities for DOT have been proposed in literature that differentiate on the kind of light modulation: continuous-wave (CW) is based on CW light, time-resolved (TR) is based on short-pulsed light, frequency-domain (FD) is based on amplitude and phase modulated illumination [8]. Among these it is possible to further add spatial modulation both on the source and detection leading to the approach of structured-light illumination and compressive sensing [9,10].

The TR approach has the important feature that time-of-flight of detected photons directly encodes space and, in a reflectance geometry where sources and detectors are placed on the same side, this means that time encodes photon penetration depth [11]. Moreover, by selecting temporal gates of the time-of-flight distribution, it is possible to directly select photons probing deeper or shallower regions of the tissue [12,13]. DOT is intrinsically an ill-posed problem and usually needs some kind of regularization that can affect the accuracy of the reconstructed optical properties [14]. This side effect leads to low-pass filtered images consequently affecting the estimate of absolute chromophore concentrations like hemoglobin, water, lipid and collagen [15]. As an example, the collagen concentration in breast lesions is fundamental for the correct estimation of breast cancer risk, and can also help lesion diagnostics [16]. The lack of accuracy is a limit for the reliability and diagnostic capability of DOT.

Most of DOT literature show reconstructed images of optical properties, where the quality reconstruction is typically evaluated as image contrast or contrast-to-noise ratio [17,18] and in terms of spatial resolution. Few works address the quantification capabilities of DOT that is the issue of accurate reconstruction of the optical properties. Works using FD-DOT in 2D have shown errors up to 15% for heterogeneities of 2.5 cm size having absorption coefficient μ_a = 0.1 cm⁻¹ over a background of μ_a = 0.05 cm⁻¹ [19]. Errors of 25% have been also reported for smaller inhomogeneities (1.7 cm size) with almost the same background and absorption perturbation [20]. In 3D FD-DOT an error of 36% has been shown for a small cylindrical object (0.8 cm diameter, 1.0 cm height) with almost the same optical properties of the above mentioned cases [21]. A more systematic work about the quantification problem in 3D FD-DOT has been reported for objects with size ranging from 1 to 2 cm and absorption perturbations up to 0.2 cm⁻¹. The error achieved, using a 3 steps reconstruction algorithm, was 27% for the 1 cm heterogeneity and 5.5% for the 2 cm one [15]. All these mentioned studies are based on FD-DOT on a cylindrical phantom with a 180 degrees arrangement of sources and detector at different planes, whereas no results are available for purely reflectance geometries, which is the most used for biomedical applications.

Despite these specific contributions, we still lack a broad consensus on shared protocols for performance assessments of DOT systems and algorithms. Such procedures and metrics would greatly help development of new instruments and reconstruction tools against objective figures, permit grading of different systems, facilitate more sound comparison of clinical results, grant quality control in clinical studies and pave the way to industrial standards. All these aspects are fundamental for a more mature growth and uptake of DOT in clinics, and to improve soundness and consistency of research.

In the more general field of diffuse optical imaging and spectroscopy, some successful attempts have been pursued in the last ten years to reach consensus in large inter-laboratory studies for performance assessment of diffuse optics instruments. We can quote, for instance, the MEDPHOT Protocol [22] for diffuse spectroscopy of homogeneous media, the BIP Protocol [23] for basic performance assessment of time-resolved systems, the NEUROPT Protocol [24] dealing with imaging in heterogeneous diffusive media. This latter is the one closest to the specificities of DOT, although it addresses more in general imaging systems and not only fully 3D tomographs.

More specifically, the NEUROPT Protocol assesses 3 key features that are sensitivity to localized small absorption changes (contrast and contrast-to-noise ratio), spatial resolution (depth selectivity and lateral resolution) and quantification of absorption changes (accuracy and linearity). In particular, accuracy and linearity are important because they are directly related to the ability of quantifying chromophore concentrations and their variations.

The idea of this paper is to make a step towards the translation of this performance assessment in TR-DOT with particular attention to the quantification problem.

In previous studies, we have already reported results on the sensitivity and spatial resolution of TR-DOT using a short source-detector separation scanning scheme [25], demonstrating that the adoption of a fast-gated single-photon avalanche diode (SPAD) has enabled the possibility to detect photons with a long time-of-flight, increasing the sensitivity at depths higher than 2.0 cm. Further, the use of a fast-gated SPAD allowed us to experimentally implement the null-distance approach with advantages in term of contrast, spatial resolution and signal throughput [26,27].

In this paper, we specifically address the issue of quantification in a TR-DOT realization in reflectance geometry. In particular, we study the effect of the source-detector separation, of the optical properties of the perturbation, and of its depth. The perturbations consist of black totally absorbing objects that – as it has been recently demonstrated [28] – are representative of a variety of perturbations with different shape and optical properties. The study is restricted to the TR approach in reflectance geometry and to purely absorbing perturbations. Phantom measurements are compared to simulations to disentangle physics or model contributions from influence of instrumentation. The reconstruction algorithm is based on the Mellin-Laplace transform (MLT) which permits to extract information in depth by time windowing the TR measurements [29].

Besides providing systematic results for the specific TR-DOT approach considered, this paper adds contributions in view of a future inter-laboratory consensus study on performance assessment of DOT instruments, which could either evolve from the NEUROPT Protocol or start as a new initiative.

The paper is organized as follows: Section 2 defines the problem, describes the experimental setup, the phantom preparation and the reconstruction algorithm. Section 3 defines the figures used to assess the quantitation performances. Section 4 displays results for linearity and accuracy both on simulation and phantom experiments. Section 5 summarizes the key findings of the study, addresses the specific factors that affect TR-DOT quantification, and discusses the implication for specific clinical problems.

2. Material and methods

2.1 Geometry

Time-domain DOT acquisitions in reflectance geometry were carried out with a horizontal scan at the surface of the phantom with a probe composed of one source and two detectors. The scan geometry was designed to obtain 30 source positions at steps of 0.75 cm with the inclusion decentered compared to the scan area to better appreciate reliability of our system to reconstruct and detect precisely in lateral (x and y) directions (Fig. 1(a)). We also investigated the influence of source detector distances by testing two probes in L configuration with one source (yellow circle) and two detectors at 1.5 cm (blue crosses) or 3.0 cm (green crosses) interfiber distances, as represented in Fig. 1(a). The center of inclusion was set at different depths z below the surface of the liquid phantom. The optodes were placed on the surface of the phantom and were held in three holes drilled in a black PVC plank. To avoid the waveguiding effect, the liquid phantom touched the black PVC holder and we removed possible air bubbles by gently dragging a finger to sweeping away bubbles. For the reconstruction, we used a mesh with a step of 0.2 cm (small grey dots in Fig. 1(a)).

 

Fig. 1 (a) Geometry of the scan x-y of both source-detector distances. The source (yellow circle) scan the yellow area. The grey disc corresponds to the position of the inclusion (x = 0 cm. y = 0 cm). The crosses are the detection positions of the couple of detector optical fibers at 1.5 cm (blue) and 3.0 cm (green) distance from the source. Each black dot is separated by 0.75 cm in both directions. (b) Instrumental setup with the phantom containing an absorbing object (detailed in section 2.2).

Download Full Size | PDF

The reconstruction is based on the analysis of the differences between the signal recorded on the inhomogeneous sample containing an inclusion and the signal recorded on a reference (homogenous medium). Such reference measurements have been acquired on an x-line scan far from the inclusion (at x = −4 cm).

2.2 Experimental setup

The experimental setup (whose schematic is reported in Fig. 1(b)) was based on a laser source providing pulses at 820 nm with a repetition rate of 40 MHz and 26 ps pulse width (Fianium Ltd, London). Light emitted from the laser was first attenuated by means of a Variable Optical Attenuator (VOA) operating from 0 up to 12 OD (Optical Density) and then injected into the medium via a 200 μm core optical fiber (NA = 0.22; 2.45 m long), as reported in Fig. 1(b).

Photons reemitted from the sample were collected in two different positions by means of two fibers (1 mm core; 0.37 NA; 2.45 m long) posed at 1.5 or 3.0 cm distance from the source, depending on the experiment. Photons harvested from each detection fiber were then focused onto a silicon fast-gated SPAD [30] (active area diameter: 100 μm) using a pair of aspheric lenses.

When a photon hits the active area of the detector, an avalanche is triggered and the fast-gated SPAD module provides a pulse that is fed as a “start” to the TCSPC board (SPC-130, Becker&Hickl GmBH, Berlin, Germany). The “stop” pulse was sent to both TCSPC boards by the laser synchronization signal. This signal was also sent to the two fast gated modules to synchronize the opening of the gate (temporal width: 5 ns). In order to enable the detection at different delays from the laser injection into the medium, the signal for the gate opening was delayed by means of a home-made programmable delayer based on transmission lines (minimum delay step: 25 ps).

The fast-gated modules were also used in the so-called “free-running” (i.e. non-gated) mode thus acquiring the full distribution of time-of-flight of re-emitted photons. In this case, the gate was opened before the first photon is reemitted and it was closed after the last photons are collected.

Both for gated and non-gated mode, at each interfiber distance, we performed the scan (following the geometry explained in the previous paragraph) using 4 different totally absorbing objects. The inclusions were posed at a depth (defined as the distance between the surface and the centroid of the inclusion) of 1.0, 1.5 and 2.0 cm. In case of gated measurements, for each scanning point and for each delay at which the SPAD is gated-ON, curves were acquired for 5 s (5 repetitions of 1 s). As required by the fast-gated acquisition technique (see Ref [31]. for details), to exploit the gating capability to collect more late photons, there is the need to increase the power injected into the phantom when increasing the gate delay. In order to guarantee a significant increase in the number of late photons when increasing the gate delay, we decided to proceed in the following way for the selection of the number of gates. We started with a first gating window opened about 2 ns before the reflectance curve peak in order to include it in the acquisition window, and we set a proper attenuation using the VOA to fit the single-photon counting statistics (i.e. the photon counting rate was kept below 1.5 Mcps, which corresponds to about 4% of the laser pulsing rate [32]) . Then, we increased the laser power injected into the phantom reducing the VOA attenuation by a factor of 5. In this situation, of course, the count rate was well above the single-photon statistical limit for TCSPC. In order to fit such limit, we then increased the delay at which the SPAD was gated-ON, thus rejecting early photons and fitting again the single-photon counting statistics. Afterwards, we repeated the procedure per each delay, up to when the maximum available power is injected into the sample. For 3.0 cm interfiber distance, we needed 3 delays while for the 1.5 cm 5 delays were necessary due to the increased number of photon reaching the detector all delays when using small source-detector separations, as explained in Ref [26]. Therefore, the acquisition times were 15 s and 25 s, respectively. In order to have the same acquisition time as for gated measurements, the collection time of photon of the non-gated acquisitions at one scanning point was 15 s and 25 s for 3.0 cm and 1.5 cm interfiber distances respectively.

2.3 Phantoms

For the realization of realistic absorption perturbations typical of biomedical situations, we followed the Equivalent Black Volume (EBV) approach, that is the use of a set of totally absorbing objects with different volumes. It was demonstrated, both with Monte Carlo simulations and phantom measurements [28], that it is possible to group different absorption perturbations of different size (volume) and intensity (absorption perturbations) in equivalence classes, whose members produce the very same effect on time-domain photon distributions over different geometries (e.g. reflectance and transmittance), source-detector distances, and background optical properties. For each class, a totally absorbing object with a given volume can be identified, yielding the same effects of all perturbations belonging to the same class. In practice, the complex combination of different possible absorption perturbations can be modelled by phantoms using a set of small black PVC cylinders with increasing volumes. Table 1 shows the dimensions and volumes (EBV) of the objects used in the present study together with their equivalent absorption perturbation (δμ_a^vol) calculated over a 1 cm³ volume sphere. These equivalent absorption perturbations were determined in [28] and correspond to the δμ_a that a spherical perturbation of 1 cm³ volume must have to produce the same perturbation of the totally absorbing object. This correspondence was validated as far as the object is not too close to the source or detector (e.g., depth < 1 cm).

Table 1. Summary of the 4 totally absorbing cylinders of various sizes with diameters, heights, corresponding volumes and the equivalent absorption perturbations in a region corresponding of a sphere of 1 cm³ [28].

View Table | View all tables in this article

The inclusions were hold in a large tank (volume 29 × 29 × 14 cm³) through a thin wire hold on the bottom of the tank and painted in white to reduce interference. The tank was filled with a water suspension of Intralipid and black ink (Higgins), yielding an absorption coefficient (μ_a) equal to 0.07 cm̵¹ and a reduced scattering coefficient (μ_s^') equal to 12 cm⁻¹ at 820 nm. We followed a standard recipe coming from the work of Spinelli et al [33].

2.4 Simulation process

To better understand the physical mechanisms of TR-DOT concerning the quantification in depth, we generated simulated measurements similar to those in the experiments (same geometric configuration). From the mathematical point of view, supposing that source and detector can be seen as points, a TR measurement is the time convolution of the Green's functions (which are intrinsic responses of the diffusive medium to a Dirac point source) and the instrumental response function (IRF) of the experimental setup. The simulation procedure involves first to convolve the experimental IRF of the SPAD with the Green's functions generated using the diffusion equation solved with our MLT approach described below. Then each curve was multiplied by a scaling factor in order to get the desired photon integral. Finally, Poisson noise comparable with experimental measurements was added. The inclusion was simulated by a 1 cm³ sphere with a δμ_a given by the Equivalence Relation reported in Table 1. Geometry, background optical properties, and reconstruction processing matched those used for the phantom measurements.

2.5 Reconstruction technique

Forward model

Light propagation in the diffusive medium was modeled using the time-domain diffusion equation:

Green’s function, noted G ($\overrightarrow{r}$,$\overrightarrow{r}$',t), is defined as the solution of Eq. (1) at position $\overrightarrow{r}$’ for a Dirac source at ${\overrightarrow{r}}_{}$. We also note G_s($\overrightarrow{r}$,t) = G(${\overrightarrow{r}}_s$,$\overrightarrow{r}$,t), G_d($\overrightarrow{r}$,t) = G(${\overrightarrow{r}}_d$,$\overrightarrow{r}$,t) and G_sd(t) = G(${\overrightarrow{r}}_s$,${\overrightarrow{r}}_d$,t) some subsets of the Green’s function where s and d are indexes for sources and detectors.

When a small perturbation on the absorption δμ_a($\overrightarrow{r}$) is applied, the Green’s functions of Eq. (1) are known to vary according to Eq. (2):

Inverse problem

Our reconstruction method is based on the work by Puszka et al [34] which showed that by combining perturbed $M_{sd}^B(t)$ and reference $M_{sd}^A(t)$measurements acquired on two configurations without (A) and with the inclusion (B) with the known Green’s functions (G^A) of configuration (A) and the estimation at iteration k of the Green’s functions (G^B(k)) of configuration (B), we obtain an equation which links measurements to the update δμ_a to be applied to the unknown absorption map (μ_a^{(k + 1)}₌ μ_a^(k)₊ δμ_a^(k)):

The advantage of such an equation is that it does not require the knowledge of the IRF of the acquisition system.

Discretization of the problem

To numerically solve the problems (forward problem, i.e. the determination of Green’s function and inverse problems, i.e. the μ_a update computation), discretizations in time and space are applied. Spatial discretization is obtained by using the finite volume method (FVM) which gives a finite partitioning in tetrahedra of the 8.0 × 6.6 × 3.6 cm³ medium with 24354 nodes. The time discretization is obtained by using the MLTs as described in [29]. It permits to transform the continuous TR signals f(t) to a few coefficients as shown in the following:

Update of the medium optical properties

After the problem discretization, Eq. (3) is transformed as a matrix equation:

The 3D map of the absorption coefficient update is obtained by minimizing the weighted least squares criterion χ² associated to Eq. (5) with a conjugated gradient method (5 iterations). The formula below of χ² is given in matrix form.

Ten iterations for all the cases (forward model and optical parameters updates) are performed to get the final reconstructed µ_a.

Constrained method

By incorporating geometric constraints in the reconstruction, it was demonstrated in a review [35] that prior information in near-infrared spectroscopy (NIRS) maximizes the accuracy in recovery the expected optical parameters. The constrained method implemented here consists in using R as a Dictionary matrix to force the absorbing inclusion on the expected position of the equivalent volume (i.e. a sphere of 1 cm³).

Because of the low resolution of DOT, one of the solutions is to do multimodal imaging and use a high resolution imaging like magnetic resonance imaging (MRI) to bring spatial information seen as a priori for DOT. Thus in this paper, we are exploring a priori spatial approach and see if the quantification can be improved with such an algorithm.

3. Measuring the quantitation

While the objectives assessment of sensitivity and spatial resolution in DOT have been addressed in many papers, and are related to the figures of contrast (C) or contrast-to-noise ratio (CNR) and to spatial localization and resolution, conversely the assessment of quantitation is less discussed.

To be in agreement with the use of the EBV approach, we evaluate the reconstructed δμ_a^vol integrated over a given volume because the maximum value δμ_a^max of the reconstructed absorption has non-physical meaning since it will depend on the effective volume of the perturbation. This approach is consistent with the findings of the EBV study. Since plenty of combinations of absorption changes and volumes yield the very same effect on the measurements, then it is not possible to assess δμ_a alone. Rather, we can estimate the equivalent δμ_a^vol corresponding to a given volume. It is not even simply the product of δμ_a and volume that is retrieved since the Equivalence Class implies a non-linear relation. Clearly, these equivalence relations are in force only for small objects (e.g. volume ≤ 1 cm³) and are ultimately related to the poor spatial resolution of DOT [28].

In this paper, we quantified the absorption variation (δμ_a) in the DOT reconstruction by comparing the integral reconstructed δμ_a,i^vol over a volume of each absorbing object i according to δμ_a,i^true which is the expected variation. In our case, the integral reconstructed δμ_a,i^vol was calculated by taking the integral over the equivalence volume of the perturbation (i.e. a 1 cm³ sphere) and by subtracting the background absorption coefficient. The background absorption coefficient was determined by taking the mean in an area without perturbation of the inclusion.

We evaluate the accuracy on the retrieval of δμ_a with a relative error ε defined as

we evaluate the linearity on the retrieval of δμ_a by fitting the dependence of δμ_a^vol vs. δμ_a^true using a 2nd order polynomial and looking at the linear and non-linear terms. More precisely, we extract the 3 parameters a, b, c using the expression:

where$y=\frac{\delta {\mu }_a^{vol}}{\gamma }$,$x=\frac{\delta {\mu }_a^{true}}{\gamma }$, and $\gamma =0.1{\text{cm}}^{\text{-1}}$ is a normalization factor to provide dimensionless coefficients. The linearity of the retrieval can be assessed considering the slope of the interpolating polynomial that is the first derivative of the previous expression:

thus, $a$ represents the non-linear distortion that must be referred to the $b$ value. The deviation from linearity can be expressed as:

Another important parameter of quantification in clinical diagnosis is precision. This parameter is typically addressed under the sensitivity framework, as contrast-to-noise ratio, which indicates how the detected perturbation stands out of noise. We have not studied this figure systematically for all combinations of δμ_a, z, ρ, and gated modality. Rather, we performed some overnight repeated scans for one inclusion and we obtained relative standard deviations less than 8% at different depths, which demonstrates a good stability of the quantified values.

4. Results

The simulation and experimental results obtained with non-gated SPAD of the reconstructed 3D absorption maps are displayed in Figs. 2 and 3 with two cut views along the transversal slice (z-y) at x = 0 and the horizontal slice (z-y) at the expected depth. The common colorbar for simulation and experiment shows the quantitative scale of the reconstructed absorption coefficient distribution.

 

Fig. 2 Numerical and experimental 3D reconstruction absorption maps with non-gated SPADs at 3.0 cm source-detector distance represented in a vertical (z-y) and a horizontal (y-x) slices both passing through the expected center of the inclusion (x = 0, y = 0 and z = 1.0 cm). From top to bottom, the 4 inclusion sizes are shown (see Table 1 for equivalent δμ_a). The black circle corresponds to a 1 cm³ sphere and is centered on the expected position of the inclusion.

Download Full Size | PDF

 

Fig. 3 Numerical and experimental 3D reconstruction absorption maps with non-gated SPADs at 3.0 cm source-detector distance represented in a vertical (z-y) and a horizontal (y-x) slices both passing through the expected center of the inclusion “size 2” (equivalent δμ_a = 0.087 cm⁻¹). From top to bottom, the 3 depths are shown. The black circle corresponds to a 1 cm³ sphere and is centered on the expected position of the inclusion (x = 0, y = 0 and given depth).

Download Full Size | PDF

The absorption perturbation appears as a spot surrounded by a quite homogenous background whose μ_a is close to the expected value of 0.07 cm⁻¹ with less than 13.1% relative error. For all the maps, a good localization in x-y and in depth z with an accuracy better than 0.2 cm is obtained. In Fig. 2, a gradual increase of the absolute absorption coefficient is observed from inclusion “size 1” (δμ_a = 0.056 cm⁻¹) to inclusion “size 4” (δμ_a = 0.37 cm⁻¹). The absorption values of the experimental results are slightly higher than in simulation. Figure 3 shows the effect of the depth on quantification. The reconstructed absorption of a given inclusion (“size 2” is shown here) decreases with the depth (1.0-1.5-2.0 cm). The experimental and simulation 3D maps are comparable and follow the same trend. In both figures (Fig. 2 and Fig. 3), artefacts appear for high absorption inclusions (white hollows below and around the absorption spot) probably due to the L configuration of the probe affecting the reconstruction. With gated SPAD, we get the same remarks on the 3D reconstruction absorption maps (not shown).

Figure 4 synthesizes the quantified absorption variations δμ_a^vol from the 3D reconstruction absorption results for all the configurations: at 1.5 or 3.0 cm interfiber distance, with non-gated or gated SPAD, for simulations and for experiments (dotted and continuous lines, respectively) at 1.0, 1.5 or 2.0 cm depth and for each inclusion. For example, the continuous and dotted brown curves in box “Non Gated” and “ρ = 3.0 cm” are absorption variations extracted from the 3D reconstruction maps of the absolute absorption distribution of Fig. 2. All the plots are increasing in accordance with the expected black curve. Concerning the linearity, the plots seem to be fairly linear up to δμ_a = 0.15 cm⁻¹ in all cases. The gated SPAD slightly improves the quantification for ρ = 1.5 cm. The shorter source-detector separation provided slightly better results in general. A possible reason for this is due to stricter photon confinement at the shorter as compared to a larger distance yielding a better contrast and thus a better quantification. We observe that the quantification decreases with depth and this effect is the same for both experiment and simulation data. This reinforces the reliability of the direct model based on the diffusion approximation. Only the experimental results with the inclusion “size 4” (δμ_a = 0.37 cm⁻¹) with ρ = 1.5 cm is far from the associated simulation which is possibly due to the limits of the conditions of EBV approach for short distances (inclusion vs surface and source-detector couple) or to the higher complexity of measurements at a short ρ.

 

Fig. 4 Overview of quantified plots of δμ_a^vol, the integral variation in absorption coefficient over a 1 cm³ spherical volume, in the different configurations: with small or large source-detector distance (ρ) (columns), by using two non-gated or gated SPADs (rows), in simulation and experiment, at 1.0, 1.5 or 2.0 cm depth and with each inclusion. The black curves correspond to the expected values δμ_a of each inclusion. The dotted curves correspond to the simulation curves and the continuous curves correspond to the experimental curves.

Download Full Size | PDF

The accuracy is obtained only at 1.0 cm up to 0.15 cm⁻¹ with relative errors inferior to 30% in the experiment (brown curves in Fig. 4). Results in simulation for deep inclusion are similar for 1.5 or 3.0 cm interfiber distance. The decrease of quantification with depth and for high absorptions may be due to the loss of resolution of the reconstruction with depth and the difficulties to get marked inclusions without smoothing of the reconstructed data. The current algorithm seems to have limits to reconstruct absorbing objects with high absorption variation above 0.2 cm⁻¹ but it is an acceptable absorption range for the target medical applications of diffuse optics imaging [16].

With the spatially constrained method (Fig. 5), the quantified plots are more linear and the reconstructed values of the absorption perturbation are much larger and closer to the expected ones for all depths. In simulation, we recover exact absorption of the background and inclusion for each depth, each mode and for both source-detector distances. With 3.0 cm interfiber distance, the experimental curves are linear with a small offset. For ρ = 1.5 cm, some problems of accuracy are visible for high absorption (δμ_a = 0.37 cm⁻¹). The constrained method which solves the depth-dependent decrease in retrieved absorption can amplify also errors or artefacts. The error of quantification could be in the high perturbation produced by the object. For short ρ this produces also a great contrast well beyond the small perturbation approximation. Thus, the error could derive by some failure of the model to reproduce correctly large perturbations. By comparing the two modes, the gated mode gives slightly better accuracy for ρ = 1.5 cm than the non-gated mode but no difference is observed for ρ = 3.0 cm.

 

Fig. 5 With constrained method, quantified plots of δμ_a in the different configurations: with small or large source-detector distance (ρ) (columns), by using two non-gated or gated SPADs (rows), in simulation and experiment, at 1.0, 1.5 or 2.0 cm depth and with each inclusion. The black curves correspond to the expected values δμ_a of each inclusion. The dotted curves correspond to the simulation curves and the continuous curves correspond to the experimental data.

Download Full Size | PDF

This method gives good perspectives for the use of TR-DOT to quantify though it still requires to know the size and position of the absorption perturbation (with another imaging modality for instance). By constraining the reconstruction on specific regions the dimension of the space of unknowns is reduced. Consequently, the problem of quantification is no more ill-posed as demonstrated by the important improvements of results in Fig. 5.

Figure 6 displays the percent of relative errors of δμ_a^vol (calculated with Eq. (7)) between the simulation and theory and between experiment and theory. The same is reported in Fig. 7 but applying the constrained method. Using the standard code (Fig. 6), relative errors are increasing with depth and with absorption perturbations independently of the gating-modality and the source-detector distance both for phantom and numerical realizations. For example, we get with the inclusion “size 2” in experiment a relative error ε = −16% on average with a standard deviation σ = 3% at 1.0 cm depth and ε = −49% on average with a standard deviation of 4% at 1.5 cm depth. For inclusion “size 3” at 1.0 cm, we obtain ε = −28% with σ = 6%. Thanks to the constrained method (Fig. 7), we recover in simulations the exact absorption variations where the average ε = 0.2% going from 0 to 4% with σ = 1%. Relative errors of phantom experiments in Fig. 7 are almost all positives and for ρ = 3.0 cm, average and standard deviation (respectively 44% and 14%) of relative errors are lower in absolute and less spreading than in Fig. 6 (respectively −53% and 22%) without the a priori approach. At the 3 depths, no significant difference is observed between non-gated and gated mode.

 

Fig. 6 Relative errors in percent (%) of δμ_a^vol between (a) simulation and theory (Sim/Th) and (b) experiment and theory (Exp/Th) calculated with Eq. (7) for each reconstruction in each configuration (ρ = 1.5 or 3.0 cm and Non Gated (NG) or Gated (G) SPAD). Colors encode depth while the filling gated or not and the color intensity the ρ.

Download Full Size | PDF

 

Fig. 7 [Constrained method] Relative errors in percent (%) of δμ_a between (a) simulation and theory (Sim/Th) and (b) experiment and theory (Exp/Th) calculated with Eq. (7) for each reconstruction in each configuration (ρ = 1.5 or 3.0 cm and Non Gated (NG) or Gated (G) SPAD). Colors encode depth while the filling Gated or Not Gated and the color intensity the ρ.

Download Full Size | PDF

To extract quantitative information on the linearity of the reconstruction, we report in Table 2 the coefficients of the polynomial fit of δμ_a^vol vs. δμ_a^true divided by a reference γ = 0.1 cm⁻¹ as defined in Section 3. The slope and the NL term are estimated for x = 1 and Δx = 1, that is around γ = 0.1 cm⁻¹. For what concerns the simulations, we observe a general trend with respect to the depth z, substantially similar for both source-detector distances and for the gated and non-gated modalities. The slope is close to 50%, 30% and 15% for z = 1.0, 1.5, and 2.0 cm respectively. In practical terms, all this means that compared to the ideal slope = 100% there is a z-dependent decrease in slope by a factor of 2, 3 and 6 for z = 1.0, 1.5, and 2.0 cm respectively. This is a strong effect, still not so huge to mask deep changes due to more superficial alterations. The NL coefficient accounts for the nonlinearity for increasing δμ_a^true and is around −20%, at all depths. Thus, the non-linearity with the increased δμ_a^true is acceptable at least for absorption changes in the order of 0.1 cm⁻¹. The c coefficient that is the offset at δμ_a^true = 0 is substantially negligible (~5%) (Table 2).

Table 2. Coefficients of the polynomial fit of reconstructed δμ_a^vol vs. δμ_a^true as defined in Section 3 for Non Gated (NG) and Gated (G) SPAD modes.

View Table | View all tables in this article

The coefficients (Table 2) related to the experiments are more scrambled, as expected, since only 4 absorption points affected by experimental noise are possibly not enough to get a robust 3-parameter fit. The general trend for ρ = 3.0 cm is substantially similar to what observed on simulations still with a higher non-linearity (more around 30%). For ρ = 1.5 cm there is more discrepancy with simulations, clearly due to noisy measurements particularly for the non-gated detector as observed in Fig. 4. Upon applying the gating, data are slightly more regular, in particular at the larger depth (z = 2.0 cm).

Table 3 displays the same parameters yet for the constrained method. On simulations the agreement is perfect with basically only b = 100% and all other terms negligible. This means that the knowledge of the exact location and size of the perturbation can completely cure the depth and absorption related decrease in the retrieved δμ_a. These results are substantially confirmed also on experiments at ρ = 3.0 cm. For ρ = 1.5 cm results are substantially altered with a strong non-linearity and non-systematic alterations. It looks like measurements at ρ = 1.5 cm are more critical, and the constrained method while fixing model-based inaccuracy, at the same time enhances any experimental artefacts. The gating does not help here.

Table 3. [Constrained method] Coefficients of the polynomial fit of reconstructed δμ_a^vol vs. δμ_a^true as defined in Section 3 for Non Gated (NG) and Gated (G) SPAD modes.

View Table | View all tables in this article

5. Discussion and conclusions

In this paper, we have addressed the quantification of DOT using time-domain reflectance measurements and Mellin-Laplace analysis both on simulations and on phantoms. The issues of sensitivity (i.e., minimum detectable focal perturbation as a function of depth) and localization (i.e., correct retrieval of lateral and depth position of the perturbation) have already been addressed elsewhere [25–27] and results reported here basically confirm those findings. Rather, this paper is focused on the quantification of the value of the reconstructed absorption perturbation.

Five main conclusions can be drawn from the results:

Taken as a whole, these results are quite encouraging since they demonstrate that for a fixed depth – e.g. in brain functional imaging at the brain cortex – absorption linearity for limited yet realistic absorption changes is preserved. This feature is important for instance in functional brain imaging or in the study of brain connectivity, since it permits to follow temporal evolutions of the signal during the exercise, or to perform spectral analysis with low distortion. The depth-decrease in the reconstructed δμ_a is clearly present and must be taken into account when comparing absolute reconstruction at different depths – for e.g. characterization of breast lesions in reflectance geometry. Still, this effect is somewhat a constant trend, not much dependent on the other parameters such as the measurement geometry, and could be somehow taken into account.

The origin of the depth- and absorption- related underestimation in δμ_a – already observed in different DOT papers – seems to be indeed due to the general spread of the reconstructed δμ_a – possibly due to the ill-posedness of the problem – leading to a dilution of the absorption change. Constraining the region of the optical perturbation – as recalled above at item #4 – completely solves this problem. In practical terms, this approach is not unrealistic, since can be implemented in co-registration with a different imaging modality yielding the size and location of the activation or suspect lesions.

Clearly, the present paper is not exhaustive since it leaves untouched the effect of other parameters, such as the background optical properties, different source-detector arrangements, and other reconstruction schemes. Still, the aim here is more on one side to appreciate the effective quantification capability of a practical DOT system, and whether this is acceptable for the specific applications, on the other to contribute to the proposal of tests and figures-of-merits to measure quantitation of δμ_a towards the proposition of shared protocols for objective performance assessment of DOT systems and algorithms.