## Abstract

Simulations and phantom measurements are used to evaluate the ability of time-domain diffuse optical tomography using Mellin-Laplace transforms to quantify the absorption perturbation of centimetric objects immersed at depth 1-2 cm in turbid media. We find that the estimated absorption coefficient varies almost linearly with the absorption change in the range of 0-0.15 cm^{−1} but is underestimated by a factor that depends on the inclusion depth (~2, 3 and 6 for depths of 1.0, 1.5 and 2.0 cm respectively). For larger absorption changes, the variation is sublinear with ~20% decrease for δμ_{a} = 0.37 cm^{−1}. By contrast, constraining the absorption change to the actual volume of the inclusion may considerably improve the accuracy and linearity of the reconstructed absorption.

© 2016 Optical Society of America

## Corrections

14 November 2016: A correction was made to the author affiliations.

## 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

^{−1}over a background of

*μ*= 0.05 cm

_{a}^{−1}[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

^{−1}. 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)).

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

^{3}volume sphere. These equivalent absorption perturbations were determined in [28] and correspond to the

*δμ*that a spherical perturbation of 1 cm

_{a}^{3}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).

The inclusions were hold in a large tank (volume 29 × 29 × 14 cm^{3}) 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̵

^{1}and a reduced scattering coefficient (

*μ*) equal to 12 cm

_{s}^{'}^{−1}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^{3} 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:

*c*is the speed of light in the medium depending of the optical index fixed at

*n = 1.4*.

*D*($\overrightarrow{r}$) is the spatial distribution of the diffusion coefficient which is defined by $1/(3{\mu}_{s}^{\text{'}}(\overrightarrow{r}))$ with${\mu}_{s}^{\text{'}}(\overrightarrow{r})$ the reduced scattering coefficient.

*μ*($\overrightarrow{r}$) is the spatial distribution of the absorption coefficients. $S(\overrightarrow{r},t)$ is the spatial and temporal distribution of the light source. $\varphi (\overrightarrow{r},t)$ is the photon density. To take into account the boundary constraints of the surface, we apply the modified Robin boundary condition. For the reconstruction, we will only take ${\mu}_{a}(\overrightarrow{r})$ as the unknown and ${\mu}_{s}^{\text{'}}(\overrightarrow{r})$ is chosen constant and equal to 12 cm

_{a}^{−1}.

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*,t) = G(${\overrightarrow{r}}_{d}$,$\overrightarrow{r}$,t) and

_{d}($\overrightarrow{r}$*G*,${\overrightarrow{r}}_{d}$,t) some subsets of the Green’s function where

_{sd}(t) = G(${\overrightarrow{r}}_{s}$*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*) of configuration (B), we obtain an equation which links measurements to the update

^{B(k)}*δμ*to be applied to the unknown absorption map (

_{a}*μ*):

_{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

^{3}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:

*n*(integer) is the order of the transform (growing from 0 to N, here N being fixed to 20) and

*p*(a positive real number) is the analysis precision set here to 3 ns

^{−1}, meaning that 3 coefficients are extracted per nanosecond. The MLT, which performs windowings on TR signal, is suitable to extract pieces of information from late photons and, therefore, to improve the quality of reconstruction of deep layers of scattering media. The ability to detect an inclusion in depth was studied in Ref [34].

### Update of the medium optical properties

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

where*W*is the sensitivity matrix and

*m*index of the volume of node.

The 3D map of the absorption coefficient update is obtained by minimizing the weighted least squares criterion *χ ^{2}* associated to Eq. (5) with a conjugated gradient method (5 iterations). The formula below of

*χ*is given in matrix form.

^{2}*L*(Left preconditioning) matrix in Eq. (6) is a diagonal matrix filled with the inverse of standard deviation of the noise on

*Y*, derived from an assumption of photon noise on the original measurements${M}_{sd}^{A}$and ${M}_{sd}^{B}$. We also introduce an

*R*(Right preconditioning) matrix to reexpress the problem on an adapted basis ([

*δμ*= RX]). In the following,

_{a}*R*is either used to attribute weight to voxels when it is a diagonal matrix. In this case, we set its elements

*R*to the square of the distance between the position of node

_{mm}*m*and the proximal source and detector locations to reinforce sensitivity in deep layers. In another case

*R*can also force identical optical parameter per predefined region.

*R*is then a “Dictionary matrix” whose columns gather voxels together into these predefined regions.

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^{3}).

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

*δμ*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}^{max}*δμ*alone. Rather, we can estimate the equivalent

_{a}*δμ*corresponding to a given volume. It is not even simply the product of

_{a}^{vol}*δμ*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

_{a}^{3}) 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

*δμ*over a volume of each absorbing object

_{a,i}^{vol}*i*according to

*δμ*which is the expected variation. In our case, the integral reconstructed

_{a,i}^{true}*δμ*was calculated by taking the integral over the equivalence volume of the perturbation (i.e. a 1 cm

_{a,i}^{vol}^{3}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

*δμ*vs.

_{a}^{vol}*δμ*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:

_{a}^{true}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:

which is the fractional deviation from a linear behavior over a spanned range of absorption of$\Delta x$. For example, referring to the numbers which will be identified in the Results session for*z*= 1 cm, a nonlinear coefficient$a=10\%$ combined with $b=100\%$ implies that a distortion from linearity (i.e. a relative decrease in slope) of $-20\%$ is expected spanning an absorption range of $\Delta x=1$ (i.e. over a range of absorption change of 0.1 cm

^{−1}). The linear coefficient $b$ should be as close as possible to 100% so as to reproduce the correct variation in

*δμ*. For high values of

_{a}*NL*, the slope must be evaluated for the effective

*x*value. A lower value of $b$ – accompanied by a low non-linearity

*NL*– indicates still a linear trend but with a reduced slope on the retrieved

*δμ*. Obviously, for $b\to 0$ and$NL=0$, the system is linear but absorption variations are so low that cannot possibly be detected. The coefficient $c$ displays the offset of the retrieval for$x=0$, thus for a null absorption perturbation. The closeness of $b$ to 100%, accompanied by a low $a$ and $c$ coefficients can be used as another measure of

_{a}*accuracy*.

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.

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

^{−1}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” (

*δμ*= 0.056 cm

_{a}^{−1}) to inclusion “size 4” (

*δμ*= 0.37 cm

_{a}^{−1}). 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

*δμ*= 0.15 cm

_{a}^{−1}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” (

*δμ*= 0.37 cm

_{a}^{−1}) 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

*ρ*.

The accuracy is obtained only at 1.0 cm up to 0.15 cm^{−1} 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^{−1} 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

^{−1}). 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.

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.

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.

*δμ*divided by a reference

_{a}^{true}*γ*= 0.1 cm

^{−1}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

^{−1}. 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

*δμ*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

_{a}^{true}^{−1}. The

*c*coefficient that is the offset at

*δμ*= 0 is substantially negligible (~5%) (Table 2).

_{a}^{true}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.

## 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:

- 1) The main parameter affecting the correct reconstruction of
*δμ*is the depth_{a}^{vol}*z*. A depth-dependence underestimation of the absorption is observed, with a reduction of the slope by a factor of 2, 3 and 6 for*z*= 1.0, 1.5 and 2.0 cm respectively. - 2) The reconstructed
*δμ*is fairly linear with respect to the increase in the real_{a}^{vol}*δμ*with the absorption change in the range 0-0.15 cm_{a}^{−1}. For higher absorption changes, a deviation from linearity is observed with a non-linear coefficient*NL*of around 20% for a change in*δμ*= 0.1 cm_{a}^{−1}(for a 1 cm^{3}volume). - 3) The adoption of a constrained approach, where the perturbation location and volume are fixed
*a priori*, completely cures depth- and absorption- reduction in the reconstructed*δμ*on simulations, and greatly improves the outcome on experiments._{a} - 4) The geometry (source-detector distance
*ρ*), as well as the adoption of a fast gating approach to suppress early photons have marginal effects both on simulations and experiments. While substantial improvements in sensitivity and localization were demonstrated adopting a short-distance [27], the use of a fast-gated approach seems to provide only a minor advantage for the issue of quantification. A possible explanation results from worse photon confinement at the larger ρ as compared to a shorter distance yielding a deterioration of contrast and thus affecting the quantification. Conversely, results at the shorter*ρ*= 1.5 cm distance are more scrambled and more prone to experimental alterations. The fast-gating does not help here. - 5) Phantom measurements substantially agree with simulations, at least for the general trends. Apart from some intrinsic higher variability, phantom measurements show a systematic small overestimation of the reconstructed
*δμ*particularly for shallower inclusions. This effects has still to be fully explained and could well reside in some experimental inaccuracies, yet, the main conclusions of items 1, 2, 3, and 4 are fully confirmed by phantom measurements._{a}

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

*δμ*– 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.

_{a}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.

## Funding

LASERLAB-EUROPE (grant agreement no. 654148, European Union's Horizon 2020 research and innovation programme).

## Acknowledgments

Portions of this work were presented at the Biomedical Optics Congress in 2016, “Quantification of effective absorption perturbations for Time-Resolved Diffuse Optical Tomography with totally absorbing objects”.

## References and links

**1. **J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys. **30**(2), 235–247 (2003). [CrossRef] [PubMed]

**2. **A. Pifferi, A. Farina, A. Torricelli, G. Quarto, R. Cubeddu, and P. Taroni, “Review: Time-domain broadband near infrared spectroscopy of the female breast: a focused review from basic principles to future perspectives,” J. Near Infrared Spectrosc. **20**(1), 223–235 (2012). [CrossRef]

**3. **Z. Yuan, Q. Zhang, E. S. Sobel, and H. Jiang, “Tomographic x-ray-guided three-dimensional diffuse optical tomography of osteoarthritis in the finger joints,” J. Biomed. Opt. **13**(4), 044006 (2008). [CrossRef] [PubMed]

**4. **J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab. **23**(8), 911–924 (2003). [CrossRef] [PubMed]

**5. **P.-Y. Lin, K. Hagan, A. Fenoglio, P. E. Grant, and M. A. Franceschini, “Reduced cerebral blood flow and oxygen metabolism in extremely preterm neonates with low-grade germinal matrix- intraventricular hemorrhage,” Sci. Rep. **6**, 25903 (2016). [CrossRef] [PubMed]

**6. **A. T. Eggebrecht, S. L. Ferradal, A. Robichaux-Viehoever, M. S. Hassanpour, H. Dehghani, A. Z. Snyder, T. Hershey, and J. P. Culver, “Mapping distributed brain function and networks with diffuse optical tomography,” Nat. Photonics **8**(6), 448–454 (2014). [CrossRef] [PubMed]

**7. **L. Di Sieno, G. Bettega, M. Berger, C. Hamou, M. Aribert, A. D. Mora, A. Puszka, H. Grateau, D. Contini, L. Hervé, J.-L. Coll, J.-M. Dinten, A. Pifferi, and A. Planat-Chrétien, “Toward noninvasive assessment of flap viability with time-resolved diffuse optical tomography: a preclinical test on rats,” J. Biomed. Opt. **21**(2), 025004 (2016). [CrossRef] [PubMed]

**8. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**(2), R41–R93 (1999). [CrossRef]

**9. **C. D’Andrea, N. Ducros, A. Bassi, S. Arridge, and G. Valentini, “Fast 3D optical reconstruction in turbid media using spatially modulated light,” Biomed. Opt. Express **1**(2), 471–481 (2010). [CrossRef] [PubMed]

**10. **J. Chen, V. Venugopal, F. Lesage, and X. Intes, “Time-resolved diffuse optical tomography with patterned-light illumination and detection,” Opt. Lett. **35**(13), 2121–2123 (2010). [CrossRef] [PubMed]

**11. **A. D. Mora, D. Contini, S. Arridge, F. Martelli, A. Tosi, G. Boso, A. Farina, T. Durduran, E. Martinenghi, A. Torricelli, and A. Pifferi, “Towards next-generation time-domain diffuse optics for extreme depth penetration and sensitivity,” Biomed. Opt. Express **6**(5), 1749–1760 (2015). [CrossRef] [PubMed]

**12. **A. Pifferi, A. Torricelli, L. Spinelli, D. Contini, R. Cubeddu, F. Martelli, G. Zaccanti, A. Tosi, A. Dalla Mora, F. Zappa, and S. Cova, “Time-Resolved Diffuse Reflectance Using Small Source-Detector Separation and Fast Single-Photon Gating,” Phys. Rev. Lett. **100**(13), 138101 (2008). [CrossRef] [PubMed]

**13. **L. D. Sieno, A. D. Mora, G. Boso, A. Tosi, A. Pifferi, R. Cubeddu, and D. Contini, “Diffuse optics using a dual window fast-gated counter,” Appl. Opt. **53**(31), 7394–7401 (2014). [CrossRef] [PubMed]

**14. **S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. **25**(12), 123010 (2009). [CrossRef]

**15. **S. Srinivasan, B. W. Pogue, H. Dehghani, S. Jiang, X. Song, and K. D. Paulsen, “Improved quantification of small objects in near-infrared diffuse optical tomography,” J. Biomed. Opt. **9**(6), 1161–1171 (2004). [CrossRef] [PubMed]

**16. **G. Quarto, L. Spinelli, A. Pifferi, A. Torricelli, R. Cubeddu, F. Abbate, N. Balestreri, S. Menna, E. Cassano, and P. Taroni, “Estimate of tissue composition in malignant and benign breast lesions by time-domain optical mammography,” Biomed. Opt. Express **5**(10), 3684–3698 (2014). [CrossRef] [PubMed]

**17. **N. Ducros, C. D’Andrea, A. Bassi, G. Valentini, and S. Arridge, “A virtual source pattern method for fluorescence tomography with structured light,” Phys. Med. Biol. **57**(12), 3811–3832 (2012). [CrossRef] [PubMed]

**18. **S. C. Davis, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Contrast-detail analysis characterizing diffuse optical fluorescence tomography image reconstruction,” J. Biomed. Opt. **10**(5), 050501 (2005). [CrossRef] [PubMed]

**19. **T. O. McBride, B. W. Pogue, E. D. Gerety, S. B. Poplack, U. L. Osterberg, and K. D. Paulsen, “Spectroscopic diffuse optical tomography for the quantitative assessment of hemoglobin concentration and oxygen saturation in breast tissue,” Appl. Opt. **38**(25), 5480–5490 (1999). [CrossRef] [PubMed]

**20. **T. O. McBride, B. W. Pogue, S. Jiang, U. L. Osterberg, K. D. Paulsen, and S. P. Poplack, “Initial studies of in vivo absorbing and scattering heterogeneity in near-infrared tomographic breast imaging,” Opt. Lett. **26**(11), 822–824 (2001). [CrossRef] [PubMed]

**21. **H. Dehghani, B. W. Pogue, J. Shudong, B. Brooksby, and K. D. Paulsen, “Three-dimensional optical tomography: resolution in small-object imaging,” Appl. Opt. **42**(16), 3117–3128 (2003). [CrossRef] [PubMed]

**22. **A. Pifferi, A. Torricelli, A. Bassi, P. Taroni, R. Cubeddu, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, J. Swartling, T. Svensson, S. Andersson-Engels, R. L. P. van Veen, H. J. C. M. Sterenborg, J.-M. Tualle, H. L. Nghiem, S. Avrillier, M. Whelan, and H. Stamm, “Performance assessment of photon migration instruments: the MEDPHOT protocol,” Appl. Opt. **44**(11), 2104–2114 (2005). [CrossRef] [PubMed]

**23. **H. Wabnitz, D. R. Taubert, M. Mazurenka, O. Steinkellner, A. Jelzow, R. Macdonald, D. Milej, P. Sawosz, M. Kacprzak, A. Liebert, R. Cooper, J. Hebden, A. Pifferi, A. Farina, I. Bargigia, D. Contini, M. Caffini, L. Zucchelli, L. Spinelli, R. Cubeddu, and A. Torricelli, “Performance assessment of time-domain optical brain imagers, part 1: basic instrumental performance protocol,” J. Biomed. Opt. **19**(8), 086010 (2014). [CrossRef] [PubMed]

**24. **H. Wabnitz, A. Jelzow, M. Mazurenka, O. Steinkellner, R. Macdonald, D. Milej, N. Żołek, M. Kacprzak, P. Sawosz, R. Maniewski, A. Liebert, S. Magazov, J. Hebden, F. Martelli, P. Di Ninni, G. Zaccanti, A. Torricelli, D. Contini, R. Re, L. Zucchelli, L. Spinelli, R. Cubeddu, and A. Pifferi, “Performance assessment of time-domain optical brain imagers, part 2: nEUROPt protocol,” J. Biomed. Opt. **19**(8), 086012 (2014). [CrossRef] [PubMed]

**25. **A. Puszka, L. Di Sieno, A. D. Mora, A. Pifferi, D. Contini, A. Planat-Chrétien, A. Koenig, G. Boso, A. Tosi, L. Hervé, and J.-M. Dinten, “Spatial resolution in depth for time-resolved diffuse optical tomography using short source-detector separations,” Biomed. Opt. Express **6**(1), 1–10 (2015). [CrossRef] [PubMed]

**26. **A. Torricelli, A. Pifferi, L. Spinelli, R. Cubeddu, F. Martelli, S. Del Bianco, and G. Zaccanti, “Time-Resolved Reflectance At Null Source-Detector Separation: Improving Contrast and Resolution In Diffuse Optical Imaging,” Phys. Rev. Lett. **95**(7), 078101 (2005). [CrossRef] [PubMed]

**27. **A. Puszka, L. Di Sieno, A. D. Mora, A. Pifferi, D. Contini, G. Boso, A. Tosi, L. Hervé, A. Planat-Chrétien, A. Koenig, and J.-M. Dinten, “Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes,” Biomed. Opt. Express **4**(8), 1351–1365 (2013). [CrossRef] [PubMed]

**28. **F. Martelli, A. Pifferi, D. Contini, L. Spinelli, A. Torricelli, H. Wabnitz, R. Macdonald, A. Sassaroli, and G. Zaccanti, “Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: Basic concepts,” J. Biomed. Opt. **18**(6), 066014 (2013). [CrossRef] [PubMed]

**29. **L. Hervé, A. Puszka, A. Planat-Chrétien, and J.-M. Dinten, “Time-domain diffuse optical tomography processing by using the Mellin-Laplace transform,” Appl. Opt. **51**(25), 5978–5988 (2012). [CrossRef] [PubMed]

**30. **G. Boso, A. Dalla Mora, A. Della Frera, and A. Tosi, “Fast-gating of single-photon avalanche diodes with 200ps transitions and 30ps timing jitter,” Sens. Actuators A Phys. **191**, 61–67 (2013). [CrossRef]

**31. **A. Tosi, A. Dalla Mora, F. Zappa, A. Gulinatti, D. Contini, A. Pifferi, L. Spinelli, A. Torricelli, and R. Cubeddu, “Fast-gated single-photon counting technique widens dynamic range and speeds up acquisition time in time-resolved measurements,” Opt. Express **19**(11), 10735–10746 (2011). [CrossRef] [PubMed]

**32. **D. O’ Connor, D.V., Phillips, *Time-Correlated Single Photon Counting* (The Royal, 1984).

**33. **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]

**34. **A. Puszka, L. Hervé, A. Planat-Chrétien, A. Koenig, J. Derouard, and J.-M. Dinten, “Time-domain reflectance diffuse optical tomography with Mellin-Laplace transform for experimental detection and depth localization of a single absorbing inclusion,” Biomed. Opt. Express **4**(4), 569–583 (2013). [CrossRef] [PubMed]

**35. **B. W. Pogue, S. C. Davis, F. Leblond, M. A. Mastanduno, H. Dehghani, and K. D. Paulsen, “Implicit and explicit prior information in near-infrared spectral imaging: accuracy, quantification and diagnostic value,” Philos Trans A **369**(1955), 4531–4557 (2011). [CrossRef] [PubMed]