Quantification in time-domain diffuse optical tomography using Mellin-Laplace transforms

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 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. 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 OCIS codes: (170.6920) Time-resolved imaging; (110.6960) Tomography; (100.3010) Image reconstruction techniques; (110.0113) Imaging through turbid media; (030.5260) Photon counting; (230.5160) Photodetectors. 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, 235–247 (2003). 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, 223–235 (2012). 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, 044006 (2008). 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, 911–24 (2003). 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 matrixintraventricular hemorrhage," Sci. Rep. (2016). 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, 448–454 (2014). 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, 25004 (2016). 8. S. R. Arridge, "Optical tomography in medical imaging," Inverse Probl. 15, R41 (1999). 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, 471–481 (2010). 10. J. Chen, V. Venugopal, F. Lesage, and X. Intes, "Time-resolved diffuse optical tomography with patterned-light illumination and detection.," Opt. Lett. 35, 2121–2123 (2010). 11. A. Dalla 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, 1749–1760 (2015). 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 SourceDetector Separation and Fast Single-Photon Gating," Phys. Rev. Lett. 100, 138101 (2008). 13. L. Di Sieno, A. Dalla Mora, G. Boso, A. Tosi, A. Pifferi, R. Cubeddu, and D. Contini, "Diffuse optics using a dual window fast-gated counter," Appl. Opt. 53, 7394–7401 (2014). 14. S. R. Arridge and J. C. Schotland, "Optical tomography: forward and inverse problems," Inverse Probl. 25, 123010 (2009). 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, 1161–71 (2004). 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, 3684–98 (2014). 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, 3811–32 (2012). 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, 050501 (2005). 19. T. O. Mcbride, B. W. Pogue, E. D. Gerety, S. B. Poplack, K. D. Paulsen, and L. O. Ulf, "Spectroscopic diffuse optical tomography for the quantitative assessment of hemoglobin concentration and oxygen saturation in breast tissue," Appl. Opt. 38, 5480 (1999). 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, 822–4 (2001). 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, 3117–3128 (2003). 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, 2104–14 (2005). 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, 86010 (2014). 24. H. Wabnitz, A. Jelzow, M. Mazurenka, O. Steinkellner, R. Macdonald, D. Milej, N. Zolek, 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, 86012 (2014). 25. A. Puszka, L. Di Sieno, A. Dalla 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–10 (2015). 26. A. Torricelli, A. Pifferi, L. Spinelli, R. Cubeddu, F. Martelli, S. Del Bianco, and G. Zaccanti, "TimeResolved Reflectance at Null Source-Detector Separation: Improving Contrast and Resolution in Diffuse Optical Imaging," Phys. Rev. Lett. 95, 078101 (2005). 27. A. Puszka, L. Di Sieno, A. Dalla Mora, A. Pifferi, D. Contini, G. Boso, A. Tosi, L. Hervé, A. PlanatChrétien, A. Koenig, and J.-M. Dinten, "Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes," Biomed. Opt. Express 4, 1351–1365 (2013). 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, 066014 (2013). 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, 5978–88 (2012). 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," Sensors Actuators A Phys. 191, 61–67 (2013). 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, 10735–46 (2011). 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, 2037–2053 (2014). 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, 569–83 (2013). 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 diagnost


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 deoxyhemoglobin 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 µ a = 0.05 cm - 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 sourcedetector 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.

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

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.nongated) 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 fastgated 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 M counts per second, 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.

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 δμ a that a spherical perturbation of 1 cm 3 volume must have to produce to obtain 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 (µ s ' ) equal to 12 cm -1 at 820 nm.We followed a standard recipe coming from the work of Spinelli et al [33].

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.

Forward model
Light propagation in the diffusive medium was modeled using the time-domain diffusion equation: ) where c is the speed of light in the medium depending of the optical index fixed at n = 1.4.D( r ! ) is the spatial distribution of the diffusion coefficient which is defined by ! is the spatial and temporal distribution of the light source.
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 as the unknown and ! ',t), is defined as the solution of Eq. ( 1) at position r !' for a Dirac source at r ! .We also note G s ( r ) 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 ( r ! ) is applied, the Green's functions of Eq. ( 1) are known to vary according to Eq. ( 2): ( ) where '*' is the time convolution product.

Inverse Problem
Our reconstruction method is based on the work by Puszka et al [34] which showed that by combining perturbed 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: The parameter 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 χ 2 is given in matrix form.
The 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 A sd M and B sd M .We also introduce an R (Right preconditioning) matrix to reexpress the problem on an adapted basis ([δµ a = RX]).In the following, R is either used to attribute weight to voxels when it is a diagonal matrix.In this case, we set its elements R mm to the square of the distance between the position of node 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.

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-tonoise 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 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 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 δµ 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 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 δµ a vol vs. δµ a true using a 2 nd order polynomial and looking at the linear and non-linear terms.More precisely, we extract the 3 parameters a, b, c using the 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 x Δ .For example, referring to the numbers which will be identified in the ).The linear coefficient b should be as close as possible to 100% so as to reproduce the correct variation in δµ a .For high values of 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 δµ a .Obviously, for 0 → b and 0 = NL , 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 , 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 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.

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" (δµ a = 0.056 cm -1 ) to inclusion "size 4" (δµ a = 0.37 cm -1 ).The absorption values of the experimental results are slightly higher than in simulation.Fig. 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).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 3 sphere and is centered on the expected position of the inclusion.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 -1 ).From top to bottom, the 3 depths are shown.The black circle corresponds to a 1 cm 3 sphere and is centered on the expected position of the inclusion (x = 0, y = 0 and given depth).
Fig. 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 -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" (δµ a = 0.37 cm -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 ρ.Fig. 4. Overview of quantified plots of δµa vol , the integral variation in absorption coefficient over a 1 cm 3 spherical volume, in the different configurations: with small or large sourcedetector 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.
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 With the spatially constrained method, 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 (Fig. 5).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 ).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. 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.Fig. 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 Table 2 (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. δµ a true divided by a reference γ = 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 δµ 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 -1 .The c coefficient that is the offset at δµ a true = 0 is substantially negligible (~5%) (Table 2).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 nonsystematic 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.

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][26][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 δµ a vol is the depth 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 δµ a vol is fairly linear with respect to the increase in the real δµ a with the absorption change in the range 0-0.15 cm -1 .For higher absorption changes, a deviation from linearity is observed with a non-linear coefficient NL of around 20% for a change in δµ a = 0.1 cm -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 δµ a on simulations, and greatly improves the outcome on experiments.
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, fast-gated approach, there seems to be only 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 δµ a 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.
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 perturbationas 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.

Fig. 1 .
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).
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 ( 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:

Results session for z = 1
distortion from linearity (i.e. a relative decrease in slope) of % 20 − is expected spanning an absorption range of 1 = Δx (i.e. over a range of absorption change of 0.1 cm -1

Fig. 6 .
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 ρ.

Fig. 7 .
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 and the color intensity the ρ.

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

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.

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.