Anne Planat-Chrétien

Michel Berger

Agathe Puszka

Laura Di Sieno

Alberto Dalla Mora

Andrea Farina

Antonio Pifferi

Tout d'abord, je souhaite remercier chaleureusement mes encadrants Jacques Derouard, mon directeur de thèse, et Lionel Hervé, mon encadrant CEA, de m'avoir guidée, suivie et accompagnée tout au long de ces trois années de doctorat. Ce fut un plaisir de travailler avec eux.

Je remercie plus spécifiquement Lionel qui m'a encadrée au quotidien en passant des heures à m'expliquer les algorithmes et en me laissant de la liberté dans l'orientation de ma thèse.

Je remercie Jean-Marc Dinten pour son

Introduction

Since the 1990's, an emerging and attractive imaging technique based on diffuse optics has been developed by researchers to probe and reconstruct in three dimensions (3D) the composition of biological tissues in depth non-invasively, non-ionizing and with potentially high specificity.

This technique is called Diffuse Optical Tomography (DOT). 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. Indeed the main constituents in biological tissues (water, lipid, melanin and hemoglobin in its oxygenated and deoxygenated forms; called chromophores) have distinct spectral signatures in the nearinfrared (NIR) range of the electromagnetic spectrum which enable to identify and quantify them. At these wavelengths, endogenous chromophores are less absorbent which allows a better penetration of light and thus the ability to probe deeper up to a few centimeters. In diffuse optics to extract information from the biological composition inside an organ, NIR light is used to illuminate the organ and the outgoing light can be collected and is then analyzed. The other clinical advantages of DOT are the flexibility and the compactness of the arrangement of the light injection and collection points. The other interests of optics in imaging diffusive medium are that light provides non-ionizing and safe radiation. Many medical applications are interested to exploit it for various organs such as breast [START_REF] Choe | Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: A case study with comparison to MRI[END_REF][START_REF] Culver | Diffuse Optical Tomography of Cerebral Blood Flow, Oxygenation, and Metabolism in Rat during Focal Ischemia[END_REF][START_REF] Enfield | Monitoring the response to neoadjuvant hormone therapy for locally advanced breast cancer using three-dimensional time-resolved optical mammography[END_REF][START_REF] Taroni | Sevenwavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification[END_REF], brain [START_REF] Eggebrecht | Mapping distributed brain function and networks with diffuse optical tomography[END_REF][START_REF] Gibson | Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate[END_REF][START_REF] Hebden | Optical tomography of the neonatal brain[END_REF][START_REF] Selb | Effect of motion artifacts and their correction on near-infrared spectroscopy oscillation data: a study in healthy subjects and stroke patients[END_REF] and others [START_REF] Khalil | Dynamic diffuse optical tomography imaging of peripheral arterial disease[END_REF] since optical properties give chemical and functional information [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF]. DOT has also the potential to characterize the vascularization of tissues by quantifying the concentration of oxy-and desoxy-hemoglobin which is useful for biomedical applications from brain activity to oncology. Recently in our laboratory, we have investigated the monitoring of the flap viability in autologous reconstruction surgery [START_REF] Di Sieno | Toward noninvasive assessment of flap viability with time-resolved diffuse optical tomography: a preclinical test on rats[END_REF] and the detection of white matter lesions in premature infant brain [START_REF] Planat-Chrétien | Time resolved optical detection for white matter lesion detection: preclinical tests on macaque brains and MRI co-registration[END_REF] into preclinical tests.

However DOT has low spatial resolution (a few millimeters) compared to conventional imaging systems because of the random photon scattering which is difficult to model. The other challenge is the penetration depth limited to a few centimeters with NIR light due to the light absorption by tissues. In optical imaging, DOT can be opposed to optical coherence tomography (OCT) which operates only at a few µm depth i.e. in superficial layers but with a high resolution of µm. This latter technology became rapidly a gold standard in ocular diseases diagnosis. The challenge of DOT is to deal with the highly diffusing nature of biological tissues in order to reconstruct in depth. Overcoming this barrier is one of the leading diffuse optical imaging research axes. Models of propagation of light have been defined in different geometries and reconstruction methods have been developed. The choice of the optimal geometry and instrumentation in accordance with clinical applications has been also studied. Another solution to tackle these drawbacks is to combine DOT with a complementary medical system; it is multimodal imaging. The main interest of combining DOT, with other clinical modalities is to gain morphological information provided by the high-resolution imaging systems in order to boost specific chemical and functional information provided by optics in the region of interest [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF][START_REF] Liang | Biomedical Optical Imaging Technologies[END_REF][START_REF] Tromberg | Assessing the future of diffuse optical imaging technologies for breast cancer management[END_REF][START_REF] Zhu | Early-Stage Invasive Breast Cancers: Potential Role of Optical Tomography with US Localization in Assisting Diagnosis[END_REF].

The focus of this PhD work is to explore the possibility to accurately characterize in 3D the composition of biological tissues in reflectance geometry using multi-wavelengths (MW) timeresolved (TR) diffuse optical tomography. Detection and spatial resolution with TR DOT were already studied in [START_REF] Puszka | Diffuse optical tomography: a time-resolved approach for reflectance measurements at short source-detector separation[END_REF]. The objective here is to move beyond simple detection, to localize, identify and quantify heterogeneities in biological tissues. The challenge is to give reliable quantitative medical images with respect to chemical and functional specificity.

In DOT, three different measurement strategies can be employed since sources and detectors can be operated into 3 different domains: continuous wave (CW) based on CW light, frequencydomain (FD) based on amplitude and phase modulated illumination and time-domain (TD) based on short-pulsed light. We choose to work in TD because of the fact that time-of-flights 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.

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 [START_REF] Di Sieno | Diffuse optics using a dual window fast-gated counter[END_REF][START_REF] Pifferi | Time-Resolved Diffuse Reflectance Using Small Source-Detector Separation and Fast Single-Photon Gating[END_REF].

The adoption of a reflectance configuration is in most cases mandatory for developing optical probes and owing that many regions of the human body are not accessible using a transmittance approach. So, the combination of the TR approach (a pulsed light source and a time-resolved detection) with a reflectance geometry enables to select the photons depending on their timeof-flights and to be specifically sensitive to superficial and deep layers in the biological tissue.

Additionally to determine the composition of biological tissues and thus separate and quantify endogenous chromophores which have specific spectra (for example: oxyhemoglobin (HbO 2 ) and desoxyhemoglobin (Hb) in Figure 1.1), multi-wavelengths acquisitions are required and the law of mixtures has to be involved in the reconstruction algorithm to retrieve the concentrations from the absorption coefficients.

To obtain endogenous chromophore features in the depth of a scattering medium, we have developed a multi-wavelengths time-resolved DOT setup combined with a 3D reconstruction algorithm. The principle of our technique is to scan a biological tissue by illuminating it with a picosecond near-infrared laser at several wavelengths and multi-positions and by collecting the backscattered light via optical fibers (at a fixed distance from the source) connected to dedicated detectors. Then reconstructions of the three-dimensional maps of chromophore concentrations of the probed tissue are obtained by using an iterative data processing technique which supposes that the propagation of light is governed by the diffusion approximation and manages Mellin-Laplace Transforms (MLT) of the MW/TR optical signals and those of a known reference medium. The quantification performances of TR DOT have been exploited with this system.

Despite a few studies addressing the issue of quantification in frequency domain [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF], the community lacks a whole protocol to perform assessments of DOT systems and especially in TD.

In this manuscript, we investigate the quantification capacity of MW TR DOT. We evaluate several parameters with might have an importance and an effect on the quantification ability of TR DOT from improvements in the reconstruction algorithm to multispectral measurements on a series of phantoms through simulation studies, a study on the shape of the instrumental response function and investigations on the quantification of the absorption coefficient in depth.

This manuscript covers the scientific work achieved during the three years of the PhD degree.

It is restricted to the TR approach in reflectance geometry of MW DOT with detection and quantification of a single object in homogeneous medium. In the perspectives section, some preliminary results on more complex media will be presented. The manuscript is organized as follows:

❼ in Chapter 1: fundamental knowledge on diffuse optics are reminded. An overview of DOT technique, the state-of-the-art of its ability to quantify and clinical transfers are introduced.

❼ in Chapter 2: we present the 3D reconstruction algorithm developed for MW TR DOT.

The whole mathematical part is explained. Different adaptations of the code (modification of the weighting of voxels (R1 code) and a priori on the geometry (a priori code or also called spatially constrained method) and data quantification analysis are defined. Simulation studies confirm that constraining the geometry improves significantly the quantification and the most useful and reliable multi-wavelengths algorithmic method is by reconstructing each wavelength measurements separately and then from the law of mixtures we obtain the concentration distribution.

❼ in Chapter 3: an evaluation to determine the limits of detection and absorption quantification of TR DOT through numerical and phantom experiments at one wavelength is detailed. Firstly, we show the impact of the shape of the instrumental response function (IRF) of the detector on the absorption contrast in depth. The decay of the IRF must be steep in order to not hide absorption information and a high dynamic range improves detectability of deep objects. Then two experimental investigations on absorption quantification, one with single-photon avalanche diodes (SPADs) and the other with new promising detectors, silicon photomultipliers (SiPMs) are reported. The experiments were carried out at the Politecnico di Milano, where these two types of detectors have been developed. To better understand the physical mechanisms of TR DOT concerning the quantification in depth, complementary numerical investigations similar to those in the experiments are given in parallel.

❼ in Chapter 4: the whole MW TR DOT algorithm is applied on a multispectral experiment campaign with a series of phantoms composed of a range of two absorbers concentrations.

We detail the choice of the composition and design of the multispectral phantoms. Material decomposition results show good gradation of the percentage of each material as expected but the accuracy decreases by a factor depending on the inclusion depth. This limit is improved by adding a spatial constraint in the reconstruction algorithm.

❼ Conclusions and perspectives: a synthesis on this PhD work and perspectives for future works are given. We also present some preliminary experimental results on bi-inclusions phantoms.

This PhD work has resulted in two publications [START_REF] Di Sieno | Time-domain diffuse optical tomography using silicon photomultipliers: feasibility study[END_REF][START_REF] Zouaoui | Quantification in time-domain diffuse optical tomography using Mellin-Laplace transforms[END_REF] in refereed journals (Biomedical Optics Express and Journal of Biomedical Optics), oral presentations and posters in international and national conferences. They are listed in Publications.

Chapter 1

Overview of Diffuse Optical Tomography (DOT)

This chapter gathers the state-of-the-art of current DOT systems and adresses the issue of quantification in scattering biological tissues. Clinical diagnostic motivations of optics in diffusive medium are highlighted, the principle of DOT is explained and a review of the endogenous quantification capacity of DOT is given. Finally, the current and future medical applications are presented and the multimodal imaging including DOT is introduced.

Optics in diffusive medium

Research on diffuse optics for medical diagnosis has been growing during the last decade thanks to the evolution of technologies such as laser, single photon counting detectors, optical fibers and computation processing speed. The medical field currently requests harmless, non-ionizing, portable, cheap and reliable medical imaging system for routine exams, screening and following cures of therapy treatments. Diffuse optical imaging has the potential to cover these challenges.

However the academic and industrial researches have some difficulties to bring out a clinical instrument with a robust software to analyze optical data. The main reasons are the problem of photon detection due to the limited penetration depth of light and the difficulty to backtrack their travels because of their random scattering in highly diffusive media.

Thanks to the improvement of the instrumentation as solid-state detector (e.g., single photon avalanche diode (SPAD)) and mathematical modelisations as Monte Carlo method, photons can now be well detected through scattering media and their travels can be estimated. To understand the limit of light penetration, we recall the physical phenomena of light-tissue interactions and the photon transport models. Then the biological interest of diffuse optical imaging for endogenous quantification in depth will be explained.

Challenging in optical penetration

Light-tissue interactions

When light meet a biological tissue, two main interactions can happen: absorption and scattering.

Light Absorption. When the laser travels through a biological medium, light can be absorbed.

The overall effect of absorption is a reduction in the intensity of the laser beam traversing the medium. This phenomenon obeys the Beer-Lambert exponential law:

I(l, λ) = I(0,λ)e -µ tot a (λ)l where I(0, λ) is the initial incident intensity, I(l, λ) is the intensity of a collimated light beam traversing a path l (cm) in a no-scattering medium and µ a is the absorption coefficient (cm -1 ) for a given wavelength λ. The absorption coefficient can thus be interpreted as the probability that a photon will be absorbed by the medium per unit length. This interaction is caused by the medium composition in endogenous chromophores -specific light absorbing molecules -.

Indeed according to the law of mixtures, the total absorption coefficient µ a,tot is the sum of the specific absorption coefficient or the molar extinction coefficient of the chromophores ε i weighted by the proportions or concentrations of the absorbers C i in presence (Eq. 1.3, i is the index of chromophore). In Figure 1.1, we can see that the spectra of oxy-and desoxy-hemoglobin (HbO 2

and Hb) which are the main chromophores in the near-infrared (NIR) spectrum are different.

In this case, optical measurements are sensitive to the oxygenation of hemoglobin. According to the wavelength, the penetration of the light is different in the tissue due to distinct spectral signatures of the chromophores and thus depends on the composition of the tissue and also on the oxygenation state. In the NIR range, biological tissues absorb less light and thus this optical window is used to probe deeper. So the relevant interest of using optics in medical applications is its high biological specificity and the possibility to quantify in depth the composition and the variation in chromophore concentration thanks to multiwavelength measurements in the NIR range.

Figure 1.1: Spectra of oxy (red line) and deoxy-hemoglobin (blue line)(HbO 2 and Hb) (Data from [START_REF] Prahl | Optical absorption of hemoglobin[END_REF]) (log scale).

Light Scattering. Scattering is caused by microscopic discontinuities in refractive index. Biological tissues have many refractive index changes. The structural properties, the different proportions and densities of lipids, proteins and other organic molecules, called in this thesis diffusive particles, are responsible of these discontinuities [START_REF] Johnsen | The Physical Basis of Transparency in Biological Tissue: Ultrastructure and the Minimization of Light Scattering[END_REF]. Theirs shape, size and refractive index perturbate the penetration of light in tissue by changing its direction. The attenuation of light due to this event is denoted by the scattering coefficient µ s which is defined as the probability of photon scattering in tissue per unit path length. Most of the time another parameter is used which combines the scattering coefficient and the anisotropy factor g (g =< cos(θ) >)

where θ is the random deviation during a scattering event) -a measure of the directionality of elastic light scattering -: it is the reduced scattering coefficient µ ′ s (µ ′ s = µ s (1 -g)). The purpose of µ ′ s is to describe the diffusion of photons in a random walk step (or transport meanfree path) of 1/µ ′ s (cm) where each step involves isotropic scattering. The value of g varies in the range from -1 to 1: g = 0 corresponds to isotropic scattering, and g = 1 to total forward scattering (Mie scattering at large particles), and g = -1 to total backward scattering. In the optical window (600-900 nm) the anisotropy factor of tissues is around 0.9 [START_REF] Cheong | A Review of the Optical Properties of Biological Tissues[END_REF][START_REF] Tuchin | Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis[END_REF].

These two coefficients (µ a and µ ′ s ) characterize the optical properties of the composition of tissues and they are wavelength dependent.

Light modelisations

Light propagation in highly diffusive media (like biological tissues) can be modeled by several means depending on the desired level of precision wanted. The most precise is the radiative tranfert equation (RTE) which is a natural description of light considered as photons. The RTE can be solved numerically by Monte Carlo (MC) stochastic method [START_REF] Wang | MCML-Monte Carlo modeling of light transport in multi-layered tissues[END_REF][START_REF] Zhu | Review of Monte Carlo modeling of light transport in tissues[END_REF], by using the full description of specific intensity L( r, θ, t) (6 dimensions, 3 for r, space, 2 for θ, 1 for t, time) or analytically with specific formula admissible for restrictive geometries [START_REF] Arridge | The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis[END_REF][START_REF] Jacques | Tutorial on diffuse light transport[END_REF]. With the MC method, photons are treated as particles, each undergoing absorption and scattering phenomena as random events. This method is commonly used but it is time consuming (it can reach weeks of calculation) since the signal-to-noise ratio (SNR) of a measurement obtained in the counting mode is:

SN R = √ N d (1.4)
where N d is the number of counted photons indexed by d the element of detection (Poisson probability distribution). Indeed, to have a good SNR and high precision it is necessary that N d be large. Large distances, complex source-detector arrangements and complex media also contribute to long time of computation. Another method, the diffusion approximation (DA) is a simplification of the RTE since it only considers photon density (ϕ( r, t)) instead of the specific intensity and it allows to have faster computations but it is only valid in some conditions: the scattering must be the dominating phenomenon (µ a ≪ µ ′ s ), all the distances must be higher than 1/µ ′ s (for example: if µ ′ s is equal to 12 cm -1 the source-detector distance has to be higher than 0.85 mm). In these conditions, the photon density ϕ( r, t) which depends on space r and time t follows the differential equation: 1 c ∂ϕ( r, t) ∂t -∇.(D( r) ∇ϕ( r, t)) + µ a ( r)ϕ( r, t) = S( r, t) (1.5) where c (c = c 0 /n) is the speed of light in the medium depending of the optical index fixed at n = 1.4 in this work. D( r) is the spatial distribution of the light diffusion coefficient which is defined by 1/(3µ ′ s ) with µ ′ s the reduced scattering coefficient. µ a ( r) is the spatial distribution of the absorption coefficients. S( r, t) is the spatial and temporal distribution of the light source.

So the photon density contains information about the optical properties of the medium (D( r) and µ a ( r)). The DA can be solved with numerical or analytical methods detailed in section 2.1.1.1.

Clinical motivations

Quantification of biological information is crucial for clinicians to give accurate clinical diagnosis. One interesting characteristic of optical imaging is its high specificity to recover biological information. For clinical applications, optical imaging has the potential to retrieve absorption or chromophore distribution inside an organ and should permit to quantify accurate concentrations of chromophores. The main chromophores in biological tissues are water (80% by volume in tissue), oxy-and desoxy-hemoglobin (HbO 2 and Hb), melanin (14.3 g/L in medium human skin) and lipid (20% in brain). This diagnostic potential relies on the fact that the absorption spectra of HbO 2 and Hb are quite different in the near-infrared wavelengths and thus it permits to do functional measurements and detect vascular activity. The absorption coefficient of the Hb is quite high in the visible (Vis) range, but is much smaller in the NIR which allows NIR light to penetrate at relatively large depth in tissues which can thus be probed. A common wavelength range called NIR therapeutic window -portion of the Vis and NIR spectrum running approximately from 600 nm to 900 nm -is adopted because biological tissues absorption is minimal; so it allows light to penetrate deeper up to several centimeters. From reference data extracted from [START_REF] Prahl | Optical absorption of hemoglobin[END_REF], Figure 1.1 shows the molar extinction coefficients of oxy-and desoxyhemoglobin as a function of wavelengths. These two endogenous chromophores have distinct where the concentrations of HbO 2 and Hb (C HbO 2 and C Hb ) in tissue can be estimated by inverting the system. From these concentrations, the oxygen saturation (SO 2 ) and the total hemoglobin concentration (C HbT ) can be obtained. These parameters are frequently used for clinical interest. They are expressed as follows:

SO 2 (%) = C HbO 2 C HbO 2 + C Hb (1.7) C HbT = C HbO 2 + C Hb (1.8)
Oxygen saturation and total hemoglobin concentration are basis parameters for pulse oximetry, monitoring and functional imaging. From normal to abnormal tissue, vascularization can change and can cause an increase of total hemoglobin concentration as angiogenesis or a decrease of oxygen saturation as tumor hypermetabolism.

One clinical challenge with DOT is to detect and quantify in depth these parameters in highly diffusive media. The objective of the thesis is part of this general aim.

In diffuse optical imaging, optical measurements are performed in the therapeutic window and light propagation is commonly modeled by the DA where the conditions are respected, because the dominating phenomenon is the scattering for NIR wavelengths; as long as the acquisition geometry respects enough distance.

The other interests of optic in imaging diffusive medium are that light provide non-ionizing and safe radiation. Technically, this technique is cheap, compact, portable and easily integrable with other imaging modalities.

N. B. : In this manuscript, we are focusing only on absorption but scattering is also an essential parameter to distinguish various biological structures and biomechanical functions. For example, collagen structure alteration related to aging skin changes the scattering properties of skin [START_REF] Calin | In vivo study of age-related changes in the optical properties of the skin[END_REF].

Methodology of DOT

To fully recover the three-dimensional (3D) chromophore distribution inside a tissue with an optical system, DOT is exploited.

DOT is a 3D imaging method to reconstruct in 3D the composition of highly diffusive biological tissues applied in depth. This modality is based on the illumination of the tissue by a set of light sources and the detection of diffusely transmitted or reflected light by a set of detectors placed on the surface of the tissue. Probing into tissues is possible up to a few centimeters. Image reconstruction to retrieve spatial and quantitative chromophore distributions is obtained by data processing the optical acquisitions. It involves a valid reconstruction model of the spatial distribution of the parameters µ a and µ ′ s from the measurements.

In the framework of this thesis, we have decided to explore the ability to quantify the endogenous chromophores of time-domain DOT using a multispectral reconstruction system employing the Mellin-Laplace transforms.

The choice of the TD approach in reflectance geometry and the reconstruction method will be explained in this section and an overview of the different DOT systems and DOT reconstruction methods will be presented. Also a non-exhaustive synthesis of the state-of-the-art of endogenous quantification in DOT will be given.

DOT systems

DOT instrumentations are composed of a collection of light sources and detectors to harvest photons and to measure the attenuation of the light intensity and/or time-of-flights of photon and of recording systems to save the acquisitions before reconstruction and analysis the data.

Depending of the clinical applications, in order to have the most adapted instrumentation the first choice is the geometry of the system (probe, helmet, etc). In diffuse optical imaging, source-detector arrangement has the advantage to be very flexible, compact and to come at 2) and ( 3) respectively. On the left, schematic distribution of photons in a sample in the reflection and transmission geometries. As a rule of thumb, the mean light penetration depth in the reflection geometry is of order ρ/2 (Extracted from [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF]).

the bedside of the patient. Two different human-machine contacts exist: the contact imaging [START_REF] Hillman | Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications[END_REF][START_REF] Taroni | Sevenwavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification[END_REF] and non-contact imaging [START_REF] Khalil | Detection of Peripheral Arterial Disease Within the Foot Using Vascular Optical Tomographic Imaging: A Clinical Pilot Study[END_REF][START_REF] Lapointe | A multi-view time-domain non-contact diffuse optical tomography scanner with dual wavelength detection for intrinsic and fluorescence small animal imaging[END_REF][START_REF] Mazurenka | Non-contact in vivo diffuse optical imaging using a timegated scanning system[END_REF]. To avoid interface coupling issues and to get high collected light efficiency, contact imaging is the most use technique. Various geometric configurations exist: planar reflection geometry, planar transmission geometry (both left Figure 1.2) and cylindrical geometry. Their uses are dependent of the medical applications. However reflectance configuration has more possible applications and is useful to design probes. Probes enable to access most anatomical sites like prostate imaging [START_REF] Laidevant | Fluorescence time-resolved imaging system embedded in an ultrasound prostate probe[END_REF], flap surgery [START_REF] Di Sieno | Toward noninvasive assessment of flap viability with time-resolved diffuse optical tomography: a preclinical test on rats[END_REF], etc.

DOT instrumentations usually differ by the type of sources: continuous-wave domain (CW), frequency-domain (FD) and time-domain (TD) light modulations.

Continuous-wave domain. CW approach uses continuous monochromatic laser source or white light to penetrate tissues. The detection can be done with one point detector like photomultipliers (PM) or photodiodes or with a camera as charge coupled device (CCD) camera. CW measurements give only information about the intensity. For one source-detector couple, the measurement is poor and very sensitive to the optical properties of superficial tissues depending of the source-detector distance. Indeed in diffuse optics, the distribution area of photons emitted from a given source and detected from a given detector is represented in Figure 1.2. It follows a shape that we commonly called "banana shape". In reflectance geometry, this representation shows a probed volume for a source-detector couple at a ρ distance.

Short source-detector distances enable to collect most of the photons which traveled in the shallow layers and had undergone few diffusion events. For larger source-detector distances, photons from superficial and deep layers are collected. In CW regime, the larger the distance between the source and the detector is, the deeper the medium can be probed but SNR is reduced. As a rule of thumb, the mean light penetration depth is of order ρ/2 in the reflection geometry. Even with a lot of of source-detector measurement couples, CW measurements do not allow µ a and µ ′ s separation and can hardly probe in depth. CW equipment is usually compact, cheap, fast and data processing is simple but suffers from poor depth sensitivity and accuracy. However some studies have developed a depth-compensation algorithm (DCA) [START_REF] Kavuri | Sparsity enhanced spatial resolution and depth localization in diffuse optical tomography[END_REF][START_REF] Niu | Development of a compensation algorithm for accurate depth localization in diffuse optical tomography[END_REF] which helps to solve the depth localization problem or providing a priori geometry information with contrast agents like fluorescence [START_REF] Zhang | Combined hemoglobin and fluorescence diffuse optical tomography for breast tumor diagnosis: a pilot study on time-domain methodology[END_REF] or with an other imaging technique like magnetic resonance imaging (MRI) [START_REF] Boas | Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function[END_REF] to constrain the reconstruction of CW DOT at the position of the heterogeneity in the medium.

Frequency-domain. FD techniques use light sinusoidally modulated in amplitude at radiofrequencies (hundreds of MHz) and record the amplitude of the output signal and the phase shift (θ) compared to the source (e.g. light-emitting diodes (LED's) with a power of approximately 100 mW) (see Figure 1.2). The phase shift brings additional information which permits to deliver information on the depth where the observed absorption changes occurred and to disentangle absorption from scattering properties of the tissue. Photodiodes or photomultipler tubes (PMT) are typically used as detectors. Compared to CW, FD brings complementary information for an equivalent acquisition time and without being expensive. We will see that in comparison with TD, the advantages of FD are that the instrumentation is cheaper and the data acquisition is faster but less information can be extracted.

Time-domain. TD approach is based on the injection of ultra-short (typically picosecond duration) laser pulses into tissue and the measurement of the pulse shape (or distribution of times of flights (DTOF) or temporal point spread function (TPSF)) (see Figure 1.2) with time-resolved (TR) detectors. There exist various kind of TR detection systems like time-gated intensified camera which consists in summing all photons arriving in a given time interval during which the detector is active and repeating this for many time intervals, streak camera or photon counting system. In our TD approach, we use photon counting system composed of TR detectors like PMs [START_REF] Hillman | Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications[END_REF][START_REF] Puszka | Time-domain reflectance diffuse optical tomography with Mellin-Laplace transform for experimental detection and depth localization of a single absorbing inclusion[END_REF], hybrid photomultipliers (HPMs), single-photon avalanche diodes (SPADs) [START_REF] Boso | Fast-gating of single-photon avalanche diodes with 200 ps transitions and 30 ps timing jitter[END_REF][START_REF] Pifferi | Time-Resolved Diffuse Reflectance Using Small Source-Detector Separation and Fast Single-Photon Gating[END_REF][START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF] or silicon photomultipliers (SiPMs) [START_REF] Mora | Fast silicon photomultiplier improves signal harvesting and reduces complexity in time-domain diffuse optics[END_REF][START_REF] Martinenghi | Spectrally Resolved Single-Photon Timing of Silicon Photomultipliers for Time-Domain Diffuse Spectroscopy[END_REF] connected to time-correlated single-photon counting (TCSPC) electronics [START_REF] Becker | The bh TCSPC Handbook[END_REF] which allows to time tag detection events. TPSF are reconstructed by doing an histogram plotting the number of photon detected in each time bin (Figure 1.3).

TD measurement includes information from a wide frequency range thanks to its time bin contrary to FD measurement which gets information from a single frequency. TPSF can be processed to give datatypes other than intensity and mean photon flight time. Equivalent CW measurements can be recovered by summing all photons independently of the time of arrival. The link between TD and FD approaches is given by the Fourier transform. In TD, the information content is potentially richer but the data acquisition is slower because the TCSPC card can registered the arrival time of one photon per source pulse (therefore, the source power must be lowered so as this condition is achieved) and the instrumentation is still expensive.

In reflectance geometry, TD regime is more appropriate to separate information coming from different layers at different depths due to the selection of photons thanks to their time-of-flights.

Photons which travel in the shallow layers are first detected because they underwent fewer diffusion events; we call them early photons. Photons which reached deep layers come back later at the surface because they have traveled a longer way; we call them late photons. So, optical properties of superficial tissue layers affect early photons and optical properties of deep layers affect late photons. An important technology advantage in TR detection is the possibility to select temporal gates of the time-of-flight distribution by deciding the gate opening of the detector at different delays and thus it increases the sensitivity in depth. By selecting successive gate openings, the full TPSF can be reconstructed with a higher dynamic range providing more sensitivity in depth. Research groups of the Politecnico di Milano were the first to develop time-gated detectors coupled to a TCSPC board for TD diffuse optics [START_REF] Di Sieno | Diffuse optics using a dual window fast-gated counter[END_REF][START_REF] Pifferi | Time-Resolved Diffuse Reflectance Using Small Source-Detector Separation and Fast Single-Photon Gating[END_REF].

At present, very few complete diffuse optics systems are available on the market for CW, for FD and especially for TD. The majority are home-made laboratory prototypes. The reasons of this restriction of the expansion of the TD technique are the current costs of the components, the intrinsic high complexity of TD diffuse optics physics and mainly the difficulty of the 3D reconstruction. However the vision of Pifferi et al [START_REF] Pifferi | New frontiers in time-domain diffuse optics, a review[END_REF] forecasts a breakthrough in the TD technology thanks to recent advancements of photonics and progress in physics of TD diffuse optics.

Reconstruction methods

Two type of methods exist in diffuse optical imaging to process the recorded data. The first one consists in exploiting the temporal data directly with NIR spectroscopy (NIRS). The second one (DOT) consists in applying a reconstruction method to recover the optical properties in reconstructed maps.

NIR spectroscopy (NIRS) method.

Some diffuse optical imaging research groups exploit directly the TPSFs to recover the composition of the medium. This approach is usually named the NIRS or functional NIRS (fNIRS).

Commonly for each source-detector couple, they recover one pixel or a collection of neighbor pixels of the whole scanning area by calculating the contrast or using a non-linear fitting procedure (as Levenberg-Marquardt algorithm or Bayesian approach) with analytical solutions [START_REF] Lapointe | A multi-view time-domain non-contact diffuse optical tomography scanner with dual wavelength detection for intrinsic and fluorescence small animal imaging[END_REF][START_REF] Liebert | Time-Resolved Multidistance Near-Infrared Spectroscopy of the Adult Head: Intracerebral and Extracerebral Absorption Changes from Moments of Distribution of Times of Flight of Photons[END_REF][START_REF] Taroni | Sevenwavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification[END_REF][START_REF] Torricelli | Time domain functional NIRS imaging for human brain mapping[END_REF] or a lookup table (LUT) inverse model based on Monte-Carlo simulations to extract optical properties. For example in [START_REF] Selb | Comparison of a layered slab and an atlas head model for Monte Carlo fitting of time-domain near-infrared spectroscopy data of the adult head[END_REF] they fitted their temporal measurements (generated with a TD NIRS system) using a Levenberg-Marquardt algorithm for nonlinear iterative least squares minimization (function lsqcurvefit from MATLAB). For this method, the medium is considered as a homogeneous semi-infinite medium or an infinite homogeneous slab [START_REF] Taroni | Sevenwavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification[END_REF] or a homogeneous two-layer slab [START_REF] Selb | Comparison of a layered slab and an atlas head model for Monte Carlo fitting of time-domain near-infrared spectroscopy data of the adult head[END_REF] for each source-detector couple and each wavelength and thus each measurement gives one pixel of the recovered map.

On the contrary, DOT consists in finding the optical properties map (unknown of the problem) which best explains the whole set of measurements. The method needs a valid model to predict measurements from the knowledge of a current optical property map and of the geometry of acquisition (this is the forward model) and a method to update the unknowns from the mismatch between predictions and actual measurements (inverse problem) (see Figure 2.1).

Generality on reconstruction method

Forward model. Usually, the diffusion approximation equation is used as forward model because in NIR DOT the conditions of validity are satisfied. By numerically or analytically solving the forward model, the measurements are predicted from the knowledge of the distribution of optical characteristics of the medium.

Inverse problem. As DOT is an indirect imaging method, the unknowns -the 3D maps of the absorption coefficients or chromophore concentrations (µ a = µ 0 a + X where µ 0 a is the absorption coefficient of a reference state and X is the variations of the optical properties between the two states) -are retrieved from measurements; it is an inverse problem. In diffuse optical imaging, the inverse problem is a nonlinear process, ill-posed and ill-conditioned. Such a problem is commonly linearized by using a perturbation approach where optical variations are considered regarding a reference state [START_REF] Arridge | Optical imaging in medicine: II. Modelling and reconstruction[END_REF]. Thus the system is written:

Y = W X (1.9)
where Y is the variations between the reference measurements and the perturbed measurements, X is the variations of the optical properties between the two states and W is the sensitivity matrix (detailed in section 2.1.2.4). Generally this system is not square and is ill-conditioned which means that the solution X is very sensitive to the noise of Y . Inverse methods such as gradient descent, conjugate gradient method, singular value decomposition and algebraic reconstruction technique can be used to solve the inverse problem by minimizing a criterion to come closer to the solution X and thus reconstruct the probed tissue.

Research on reconstruction methods for DOT started in the 1990's. Different algorithms have been developed depending on the choice of the direct model, the DOT systems, the techniques of regularization, the methods of discretization, etc. Some reviews concerning reconstruction methods in DOT detail more precisely the algorithms [START_REF] Arridge | Optical tomography: forward and inverse problems[END_REF]2].

Existing reconstruction algorithms.

Depending on the DOT system, different reconstruction methods exist.

Public software packages. Some DOT open-source software packages can be downloaded for free online for the CW, FD and TD configurations.

❼ NIRFAST (Near Infrared Fluorescence and Spectral Tomography) [START_REF] Dehghani | Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction[END_REF][START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF] ware is associated to the MONSTIR system [START_REF] Hillman | Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications[END_REF]. TOAST employs the FEM to model the propagation of light through tissue using the diffusion approximation to the radiative transfer equation. Image reconstruction involves iteratively adjusting the optical properties assigned to the FEM mesh to optimize the matching of the model to the data [START_REF] Arridge | Image reconstruction in optical tomography[END_REF][START_REF] Schweiger | The finite element method for the propagation of light in scattering media: Boundary and source conditions[END_REF].

Active

Advantages of the TD approach

After this overview of the DOT methodology, we discuss our choice of the TD approach. The advantages of the TD approach already mentioned are the rich-information from one measurement and the interesting combination with the reflectance geometry which links the layer depth with the time of arrival of photons at the detector. Indeed, for us, the best geometry to carry out DOT diagnosis is the reflectance geometry because of the many possible applications and the compactness of associated probes. One other specific advantage of TD DOT is the possibility to use different datatypes to process the TR data. The principle is to time filter (i.e. timewindowing) TR measurements for extracting information from them at different time windows.

In reflection geometry configuration, by employing time-windowing or moments analysis, it is possible to separate contribution of photons according to their arrival times and thus separates the information depending on the depth. Different datatypes (moments i.e. Mellin transform, central moments, Laplace transform, Mellin-Laplace transform) have been defined in the literature, applied and compared [START_REF] Arridge | Image reconstruction in optical tomography[END_REF][START_REF] Hervé | Time-domain diffuse optical tomography processing by using the Mellin-Laplace transform[END_REF][START_REF] Liebert | Determination of absorption changes from moments of distributions of times of flight of photons: optimization of measurement conditions for a two-layered tissue model[END_REF][START_REF] Puszka | Diffuse optical tomography: a time-resolved approach for reflectance measurements at short source-detector separation[END_REF]. A study on the comparison of different datatypes (Moments and Mellin-Laplace transforms) was carried out and was described in section 2.1.3 and 2.1.4 of Dr. A. Puszka's PhD manuscript [START_REF] Puszka | Diffuse optical tomography: a time-resolved approach for reflectance measurements at short source-detector separation[END_REF]. The conclusions of this study were that the moments and Mellin-Laplace transforms (MLT) are both interesting for TR DOT in reflectance geometry but the latter having the advantage of acting as a time-filter which enables time-windowing of the signal by reducing the time dimension and increasing the sensitivity for late time.

Review of endogenous quantification in DOT

Until recently, the aim of the assessment of DOT technology was detection, localization, sensibility and spatial resolution [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF][START_REF] Puszka | Spatial resolution in depth for timeresolved diffuse optical tomography using short source-detector separations[END_REF]. Most of the studies using DOT systems for clinical applications show reconstructed images with contrast or contrast-to-noise ratio [START_REF] Taroni | Sevenwavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification[END_REF] which aim to qualitative reconstructions and not quantitative reconstructions. This qualitative approach excludes the potential ability of DOT to quantify the tissue content although it is an important issue to give an accurate diagnosis. The reason of this lack of studies on quantification is probably due to the limited depth resolution that causes significant partial volume error and prevents absolute amplitude accuracy in the estimates of chromophore concentration [START_REF] Boas | Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy[END_REF]. Few works address the quantification capabilities of DOT, that is the issue of accurate reconstruction of the optical properties [START_REF] Boas | Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy[END_REF][START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF][START_REF] Pogue | Implicit and explicit prior information in near-infrared spectral imaging: accuracy, quantification and diagnostic value[END_REF] In literature, a large proportion of papers interested in the problem of absorption quantification is coming from researches on NIRS technique which contrary to DOT doesn't go through a reconstruction approach as explained in section 1.2.2.

Works using TD fNIRS on several standard homogeneous solid phantoms with a range of absorption coefficients (from 0.07 cm -1 to 0.28 cm -1 with a step of 0.07 cm -1 ) and a range of the reduced scattering coefficient (5, 10, and 15 cm -1 ) show retrieved values with good linearity and good accuracy with an average relative error of about <4% for the absorption and <8% for the reduced scattering coefficients [START_REF] Re | Multi-channel medical device for time domain functional near infrared spectroscopy based on wavelength space multiplexing[END_REF]. These measurements were performed in a reflectance geometry with one source and eight detectors positions on a line at a source-detector distance of 30 mm. For deep quantification, another paper [START_REF] Selb | Comparison of a layered slab and an atlas head model for Monte Carlo fitting of time-domain near-infrared spectroscopy data of the adult head[END_REF] have conducted realistic TD-NIRS measurements on a human brain atlas model with a fixed extracerebral absorption (set to 0.1 cm -1 ) and a range of brain absorption varying between 0.05 and 0.3 cm -1 by steps of 0.05 cm -1 .

The retrieved absorption of the brain is fairly accurate with a small underestimation for values larger than 0.2 cm -1 . These two papers assessed their TD NIRS system with a full study of quantification on a range of optical properties. Solutions to improve the absolute quantification in NIRS have been presented and demonstrated in a review [START_REF] Pogue | Implicit and explicit prior information in near-infrared spectral imaging: accuracy, quantification and diagnostic value[END_REF] both on tissue-simulating phantom studies and endogenous in vivo breast imaging. It concluded that spatial and spectral prior informations are important to optimize the quantitative accuracy.

Concerning DOT technique, research on this topic is limited. Most of the DOT studies show reconstruction results usually by normalizing the image scale or relative values for dynamic DOT [START_REF] Carp | Hemodynamic signature of breast cancer under fractional mammographic compression using a dynamic diffuse optical tomography system[END_REF][START_REF] Khalil | Dynamic diffuse optical tomography imaging of peripheral arterial disease[END_REF]. Also the majority of the papers carried out measurements with only one kind of absorbing object and changed other parameters like the depth [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF] for studying detection, the distance between two objects [START_REF] Puszka | Spatial resolution in depth for timeresolved diffuse optical tomography using short source-detector separations[END_REF] for estimating the spatial resolution, etc. The quantitative aspect is not systematically studied [START_REF] Corlu | Diffuse optical tomography with spectral constraints and wavelength optimization[END_REF]. In FD, one study fully evaluates the quantification performance of DOT system with a gradual range of optical variations [START_REF] Srinivasan | Improved quantification of small objects in near-infrared diffuse optical tomography[END_REF] but to our knowledge no study in TD is dedicated to this issue.

Below is a non-exhaustive list of the existing DOT studies which have investigated the problem of accurate endogenous quantification. Works using FD DOT in 2D have shown errors up to 15%

for heterogeneities of 25 mm size having absorption coefficient µ a = 0.1 cm -1 over a background of µ a = 0.05 cm -1 [START_REF] Mcbride | Spectroscopic diffuse optical tomography for the quantitative assessment of hemoglobin concentration and oxygen saturation in breast tissue[END_REF]. Errors of 25% have also been reported for smaller inhomogeneities (17 mm size) with almost the same background and absorption perturbation [START_REF] Mcbride | Initial studies of in vivo absorbing and scattering heterogeneity in near-infrared tomographic breast imaging[END_REF]. In 3D FD DOT, an error of 36% has been shown for a small cylindrical object (8 mm diameter, 10 mm height) with almost the same optical properties of the above mentioned cases [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF]. However by adding a priori spatial information, the absorption coefficient is accurately retrieved with an error reduced to 3% [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF]. A clinical 3D reconstruction of Hb and HbO 2 from the breast of a human volunteer shows not so good quantitative values but coherence increases in the tumor region compared with the background [START_REF] Dehghani | Multiwavelength threedimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results[END_REF]. A more systematic work about the quantification problem in 3D FD DOT has been reported for simulated objects with size ranging from 10 to 20 mm and absorption perturbations up to 0.2 cm -1 positioned at 10 mm from the edge. The error achieved, using a 3 steps reconstruction algorithm, was 27% for the 10 mm heterogeneity and 5.5% for the 20 mm one [START_REF] Srinivasan | Improved quantification of small objects in near-infrared diffuse optical tomography[END_REF] (see Figure 1.4(b)). All of these studies are based on FD DOT on a cylindrical phantom with a 180 degrees arrangement of sources and detector at different planes. In [START_REF] Dehghani | Multiwavelength threedimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results[END_REF], the authors have proposed modifications of the algorithm with moresophisticated regularization and addition of a priori information to improve the quantification of the reconstruction. The improvement of the quantification by forcing the reconstruction on the a priori position of the object was confirmed in [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF] and [START_REF] Srinivasan | Improved quantification of small objects in near-infrared diffuse optical tomography[END_REF] where the errors for objects of 10 mm or 20 mm diameter are less than 6% (see Figure 1.4).

About TR DOT, works in 2D on a cylindrical phantom show an error of 16% for the 5 mm diameter rod (µ a = 0.1 cm -1 ) and 8% for the background (µ a = 0.05 cm -1 ) in absorption [START_REF] Nouizi | Improvement of absorption and scattering discrimination by selection of sensitive points on temporal profile in diffuse optical tomography[END_REF]. In reflectance geometry through a solid two-layered slab phantom (µ a = 0.36 cm -1 and µ ′ s = 23 cm -1 at 780 nm) interesting semi-3D TR DOT reconstructions of a 20 mm diameter target (µ a = 0.53 cm -1 , 12 mm depth) probed at three wavelengths are presented in [START_REF] Gao | Optical tomographic mapping of cerebral haemodynamics by means of time-domain detection: methodology and phantom validation[END_REF] with a good localization and an accurate quantification with an error of 7%.

In time domain, the majority of the papers on diffuse optics imaging shows reconstruction images with normalized scale [START_REF] Taroni | Sevenwavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification[END_REF] or percentage change compared to the background [START_REF] Cooper | MONSTIR II: A 32-channel, multispectral, time-resolved optical tomography system for neonatal brain imaging[END_REF] i.e. relative quantification except the specific contributions mentioned above. To our knowledge, no results are available for reflectance geometry in TR DOT with a complete study on a gradual range of embedded objects. In this manuscript, we specifically address the issue of quantification in TR DOT in reflectance geometry and test the impact of spatial a priori information as applied in [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF] in case of FD.

Medical applications and multimodality

DOT has specific advantages such as noninvasiveness, nondestructiveness, richness of physiological information content, flexibility, compactness, which enable a wide field of medical applications. Cancer diagnosis, treatment monitoring, screening, surgery guide, functional mapping on organs are thus possible with this medical imaging device [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF][START_REF] Hillman | Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications[END_REF]. However, the low resolution limits the anatomical information which is essential for diagnosticians. To tackle this drawback, complementary medical systems can be coupled with DOT; it is multimodal imaging. This approach is in full expansion and many researches are conducted where the potential of DOT is well exploited.

Clinical applications

DOT technology has been developed for more than two decades but the translation of optical tomographic imaging methods into clinical practice is still limited and has difficulty to establish itself in the demanding healthcare market. Main systems are built in laboratory and are custommade prototypes. The industrial investment is currently limited.

A non exhaustive list of clinical applications are presented below:

Breast. Breast cancer is the leading type of cancer in women thus it is a major health issue.

The gold standard diagnosis system for breast cancer prevention is X-ray mammography which exists for 50 years. This technique is ionizing and painful for women because of the need of breast compress to reduce the dose of radiation. Mammography has a sensitivity of about 84% [START_REF]Performance Measures by Age[END_REF], so breast malignant or begnin tumors are detectable but some false positive results can still happen. DOT or NIRS could be an alternative technique with more physiological specificity and probably less discomfort for women. Up to now, the majority of the mammography system are laboratory prototypes and interesting clinical studied have been investigated by several research groups. Since 2004, some clinical NIRS studies are led by the Politecnico di Milano group [START_REF] Pifferi | Fourwavelength time-resolved optical mammography in the 680 980-nm range[END_REF][START_REF] Quarto | Estimate of tissue composition in malignant and benign breast lesions by time-domain optical mammography[END_REF][START_REF] Taroni | Sevenwavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification[END_REF][START_REF] Taroni | Time-resolved optical mammography between 637 and 985 nm: clinical study on the detection and identification of breast lesions[END_REF] on breast cancer with their MW TR optical mammograph designed to collect projection images of the lightly compressed breast in the same geometry as conventional X-ray mammography and operating beyond 900 nm including the absorption peaks of lipids and water. In [START_REF] Quarto | Estimate of tissue composition in malignant and benign breast lesions by time-domain optical mammography[END_REF], 62 patients with breast lesions were scanned in order to estimate the tissue composition of malignant and benign breast lesions by time-domain optical mammography.

The conclusions show a statistically significant discrimination where tumors have higher oxyhemoglobin and collagen content than benign lesions (Figure 1.5(a)). In another clinical study of the Politecnico di Milano [START_REF] Taroni | Is Collagen an Independent Risk Factor for Breast Cancer?[END_REF], it was pointed out that collagen content could be a dependent breast cancer risk factor with respect to the mammographic density. Another group, the Physikalisch-Technische Bundesanstalt group in Berlin has developed a mammographic system based on the same technology by employing mechanical scans of the optodes in transmittance geometry but with more restricted wavelengths [START_REF] Grosenick | Concentration and oxygen saturation of haemoglobin of 50 breast tumours determined by time-domain optical mammography[END_REF][START_REF] Grosenick | Development of a time-domain optical mammograph and first in vivo applications[END_REF]. Another approach of optical mammography with a conical geometry of the source and detector positions was undertaken. With this geometry, a FD DOT system with a reconstruction process also enables the detection of tumors and an apparent increase of oxy-deoxy hemoglobin concentration can be quantified [START_REF] Dehghani | Multiwavelength threedimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results[END_REF]. Going back to TR approach, a study on uncompressed breast with a hemispherical DOT system and a liquid-coupled interface show the possibility to scan breast without harming of patient, and the study on 38 probed women has permitted to highlight a higher blood volume in the tumor than in the surrounding tissue presumably as a result of tumor hypervascularization [START_REF] Enfield | Three-dimensional timeresolved optical mammography of the uncompressed breast[END_REF]. Another interest of DOT for breast applications is the monitoring of neoadjuvant chemotherapy (hormone therapy) [START_REF] Choe | Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: A case study with comparison to MRI[END_REF][START_REF] Enfield | Monitoring the response to neoadjuvant hormone therapy for locally advanced breast cancer using three-dimensional time-resolved optical mammography[END_REF]. Choe et al [START_REF] Choe | Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: A case study with comparison to MRI[END_REF] presents a study of the monitoring of one breast during neoadjuvant chemotherapy where 3D CW DOT results in transmission geometry are compared with MRI images. The retrieved position of the tumor correlated with the MRI images and the tumor volume and the total hemoglobin concentration contrast quantified with the dynamic DOT images decrease over the course of chemotherapy as expected. Also in TD, three case studies on monitoring patients confirmed that DOT can follow physiological changes within the tumor in response to treatment and may provide potentially valuable physiological information [START_REF] Enfield | Monitoring the response to neoadjuvant hormone therapy for locally advanced breast cancer using three-dimensional time-resolved optical mammography[END_REF].

About the oxygen saturation in the breast tumor, there exists contradictory interpretations (increase for [START_REF] Dehghani | Multiwavelength threedimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results[END_REF], decreased for [START_REF] Enfield | Three-dimensional timeresolved optical mammography of the uncompressed breast[END_REF] or no difference [START_REF] Choe | Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: A case study with comparison to MRI[END_REF]) which are probably caused by the stage or type of the tumor or the sensitivity of the system [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF]. Still some progress have to be done on this application to integrate optical mammograph in hospital. The lack of structural resolution leads some researchs to enhance contrast and spatial resolution by involving contrast agent like fluorescence [START_REF] Zhang | Combined hemoglobin and fluorescence diffuse optical tomography for breast tumor diagnosis: a pilot study on time-domain methodology[END_REF] and thus add a priori information in the reconstruction process or by incorporating other imaging modalities. The multimodal approach can potentially improve breast diagnosis and raise the quantification potential of optical mammography [START_REF] Tromberg | Assessing the future of diffuse optical imaging technologies for breast cancer management[END_REF].

On the industrial side, we identified two commercial systems [START_REF] Van De Ven | Optical imaging of the breast[END_REF]. Brain. Another important field is brain imaging: in this case a map of oxy-and deoxyhemoglobin concentrations is fundamental for the diagnosis of injuries like ischemia [START_REF] Culver | Diffuse Optical Tomography of Cerebral Blood Flow, Oxygenation, and Metabolism in Rat during Focal Ischemia[END_REF], hemorrhage [START_REF] Hebden | Optical tomography of the neonatal brain[END_REF], white matter lesion of premature infant [START_REF] Planat-Chrétien | Time resolved optical detection for white matter lesion detection: preclinical tests on macaque brains and MRI co-registration[END_REF] and for functional imaging during a variety of tasks [START_REF] Eggebrecht | Mapping distributed brain function and networks with diffuse optical tomography[END_REF][START_REF] Torricelli | Time domain functional NIRS imaging for human brain mapping[END_REF]. Brain imaging often uses NIRS or functional NIRS (fNIRS) and less DOT technique.

Currently, the most attractive brain application for research and industry is brain activity with the possibility of a whole-head mapping and compatible with implanted metal and electronic devices. Mapping hierarchical language areas with a 96-source and 92-detector high-density DOT system in comparison with functional MRI were presented in [START_REF] Eggebrecht | Mapping distributed brain function and networks with diffuse optical tomography[END_REF] and showed good overlap (Figure 1.6). A review foresees a flourishing future for TD fNIRS in the neuroimaging community [START_REF] Torricelli | Time domain functional NIRS imaging for human brain mapping[END_REF] where recent advances in photonic technologies might allow TD fNIRS to better recover biological information in deep layers of tissues and potentially to overtake CW fNIRS [START_REF] Boas | Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy[END_REF]. Detection of cerebral hemodynamics response to a breath hold task [START_REF] Selb | Functional brain imaging with a supercontinuum time-domain NIRS system[END_REF] or finger tapping task which decoded changes in deoxy-hemoglobin concentration in some areas in the brain is now possible with TD fNIRS.

Another important target is the neonatal brain imaging where some regions of skull are soft because of the fontanel of infant, which enable better penetration of light. Thus DOT has an interesting imaging profile for newborn brain imaging: it is painless, safe, non-invasive, portable, compact and light can travel in depth. In 1999, a group of researchers at University College London developed their 32-channel time-resolved optical tomographic imaging instrument prototype (so-called MONSTIR) in the aim to image the head of neonate [START_REF] Schmidt | Development of a time-resolved optical tomography system for neonatal brain imaging[END_REF][START_REF] Schmidt | A 32-channel time-resolved instrument for medical optical tomography[END_REF] and a new generation (MONSTIR II) with multispectral performances is described in [START_REF] Cooper | MONSTIR II: A 32-channel, multispectral, time-resolved optical tomography system for neonatal brain imaging[END_REF]. The device was tested on a clinical study with six infants and the first 3D reconstructed images show change in total hemoglobin during passive movement of each arm [START_REF] Gibson | Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate[END_REF]. However the low resolution (< 1 cm) affects precision in the localization of the functional information. To solve this problem, they suggested to incorporate prior anatomical information from generic MRI scans of newborn infants [START_REF] Hebden | Optical tomography of the neonatal brain[END_REF]. Hemorrhage in neonatal infant can be quantified with TR DOT (see Figure 1.5(b)) [START_REF] Hebden | Optical tomography of the neonatal brain[END_REF]. Our lab has recently investigated the detection of white matter lesions of the premature with a preclinical test on young macaque brains [START_REF] Planat-Chrétien | Time resolved optical detection for white matter lesion detection: preclinical tests on macaque brains and MRI co-registration[END_REF].

Other public health issues and questions could potentially be studied with functional DOT as stroke [START_REF] Selb | Effect of motion artifacts and their correction on near-infrared spectroscopy oscillation data: a study in healthy subjects and stroke patients[END_REF], sleep apnea, etc. The neuroimaging is an application in full expansion and a few companies are interesting and starting to exploit optical imaging. At present, functional optical imaging are mostly systems for research and they are not approved for any particular medical indication. They are in similar situation as fMRI. The majority of the commercial devices is based on NIRS technology combined with a helmet comporting lot of optical fibers for illumination and detection. The aim of these systems is the investigation of cerebral hemodynamics and the revelation of the secrets of the human brain. Some companies have developed compact, mobile and user-friendly systems [START_REF] Musa | Computational Optical Biomedical Spectroscopy and Imaging[END_REF]. In CW domain, there exist ETG 4000 from Hitachi (Japan), NIRSscout from NIRx (USA)

and OxyMon from Artiris (The Netherlands). And in FD, there is Imagent from ISS (USA). For CW DOT technique, Cephalogics worked on an equipment (still tested) with clinical interest, promising valuable measures to clinicians of cerebral perfusion status in brain-injured patients.

Vascular imaging. The strong point of DOT is the capacity to detect variation in vascularization into tissues. DOT can be a useful tool for diagnosis cardiovascular diseases or in-surgery monitoring the vascularization. Some researches groups are looking into peripheral arterial disease (PAD) which is the narrowing of arteries due to plaque accumulation and in the vascular walls. PAD leads to foot ulcers, open wounds, and gangrene when untreated. Anatomical information on the position, size and optical information can be recovered in foot thank to optical tomographic reconstruction [START_REF] Khalil | Dynamic diffuse optical tomography imaging of peripheral arterial disease[END_REF][START_REF] Zhao | Time-resolved diffuse optical tomographic imaging for the provision of both anatomical and functional information about biological tissue[END_REF] (see Figure 1.5(c)). Works on CW dynamic DOT on three clinical cases proved the capability to image the temporal response of the vasculature to an external provocation and uses this response as a diagnostic marker of PAD in foot [START_REF] Khalil | Dynamic diffuse optical tomography imaging of peripheral arterial disease[END_REF]. This research group from Columbia University has realized a prototype called VOTI (Vascular optical tomographic imaging) [START_REF] Khalil | Detection of Peripheral Arterial Disease Within the Foot Using Vascular Optical Tomographic Imaging: A Clinical Pilot Study[END_REF] which has the potential for assessing vascular disease in diabetic patients.

Recently, our lab proposed in collaboration with the Politecnico di Milano the monitoring of flap viability in autologous reconstruction surgery using TR DOT. The vascularization of these tissues is fundamental and a fast postoperative control is important for the success of the operation. A preclinical test confirms the interest [START_REF] Di Sieno | Toward noninvasive assessment of flap viability with time-resolved diffuse optical tomography: a preclinical test on rats[END_REF]. The forearm is a convenient organ where DOT can also be exploited as a diagnostic tool and some groups show good results on physiological changes in response to exercise like with TR DOT in cylindrical geometry [START_REF] Zhao | Time-resolved diffuse optical tomographic imaging for the provision of both anatomical and functional information about biological tissue[END_REF] and combined with MONSTIR system [START_REF] Hillman | Time resolved optical tomography of the human forearm[END_REF][START_REF] Hillman | Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications[END_REF].

Other applications. Diffuse optical imaging offers a broad field of applications also on other organs. In dermatology for the screening of skin cancer, the non-invasive, non-destructive and non-ionising characteristics of this technique are attractive for the imaging of the dermis and epidermis for improved differentiation between malignant and benign melanoma. A device is already commercialized by MELA Sciences and FDA approved for diagnosis of melanoma called MelaFind. It is a non-invasive tool that provides additional information to dermatologists during melanoma skin examinations. Other medical interest, the detection of joint pathologies as rheumatoid arthritis where research has been done on image classification of FD DOT images of finger joints between healthy and infected subjects [START_REF] Montejo | Computer-aided diagnosis of rheumatoid arthritis with optical tomography, Part 2: image classification[END_REF].

Concerning the industrial field, we identify a few companies for human diagnosis but they have some difficulty to obtain health and marketing authorizations. So, some companies develop optical computing tomographs for only small animal and preclinical applications for research such as Perkin Elmer, Bruker in vivo Xtreme, TriFoil Imaging.

Multimodal imaging

Recent advances in clinical imaging heads towards multimodal imaging. Multimodality imaging with two or more imaging modalities allows integration of the strengths of individual modalities, while overcoming their limitations and permits to help clinicians to have a more accurate The main interest of combining diffuse optical imaging, particularly DOT, with other clinical modalities is to gain morphological information provided by the high-resolution imaging systems in order to boost specific chemical and functional information provided by optics. The structural map acquired using the complementary modalities is used as a priori constrained for optical reconstruction process and enables great improvements in DOT reconstruction quality and quantitative biological information. This allows to overcome the limits of DOT which suffers from low-resolution and associated contrast dilution [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF][START_REF] Liang | Biomedical Optical Imaging Technologies[END_REF][START_REF] Pifferi | New frontiers in time-domain diffuse optics, a review[END_REF] .

MRI/DOT. MRI is a medical imaging technique based on nuclear magnetic resonance of atoms within the body induced by the application of radio waves. It provides high-resolution images of the anatomy of the body using strong magnetic fields, radio waves, and field gradients. The combination of MRI/DOT is interesting because structural MRI images can guide the DOT reconstructions. Segmentation of different tissue types and 3D mesh creation from MRI images can be processed by a software tool such as NIRFAST. Then this spatial prior knowledge is injected into the DOT inverse problem in order to help the solving and it significantly improves the retrieval of the optical properties [START_REF] Boas | Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function[END_REF][START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF]. The MRI/DOT fusion is technologically challenging because of the compatibility with MRI constraints. The DOT instrumentation is often composed of long optical fibers to operate outside the MRI room to avoid interference of the strong magnetic field with the electronic part of DOT system [START_REF] Liang | Biomedical Optical Imaging Technologies[END_REF].

CT/DOT. CT is based on X-ray imaging technique in which cross-sectional images are reconstructed by computer-processed combinations of many X-ray images taken from different angles. For the combination CT/DOT, breast imaging using X-ray mammography is an interesting clinical application [START_REF] Tromberg | Assessing the future of diffuse optical imaging technologies for breast cancer management[END_REF]. The instrumentation design is also difficult with CT to coregister tomographic x-ray and optical imaging which limit the clinical translation [START_REF] Liang | Biomedical Optical Imaging Technologies[END_REF]. A research group at Massachusetts General Hospital realized a hybrid optical and X-ray mammography imaging system which integrated an optical system with a digital breast tomosynthesis.

Initial clinical results led to a better understanding of the relationship between the functional contrast provided by optical imaging and the structural contrast afforded by X-ray imaging [START_REF] Zhang | Coregistered tomographic x-ray and optical breast imaging: initial results[END_REF].

US/DOT. US is a contact and non-ionizing system composed of a handheld probe which transmits and receives high frequency acoustic waves in the order of several megahertz. Ultrasound images are formed by the transmission of ultrasound waves, the detection of echoes from the layers between different tissues and the interpretation of these echoes. The merging of the two instruments US and DOT is easier because of the flexibility and compactness of the ultrasound probe. Combining US/DOT leads to restricted geometry in a reflectance geometry approach [START_REF] Liang | Biomedical Optical Imaging Technologies[END_REF]. Clinical trial on early-stage breast cancers was investigated on a cohort of women and results demonstrated that small indeterminate breast lesions localized with US alone can be accurate diagnosis (between malignant and benign lesions) thanks to quantitative information of the total hemoglobin concentration provided by NIR technique [START_REF] Zhu | Early-Stage Invasive Breast Cancers: Potential Role of Optical Tomography with US Localization in Assisting Diagnosis[END_REF]. Figure 1.7 shows the coregistered US/NIR system and some clinical results of a benign and a malignant breast lesions.

PAT/DOT. PAT is a novel optical medical imaging composed of a light delivery system and an ultrasound receive system. This technique is based on photoacoustic effect which is the formation of acoustic waves following light absorption in biological tissues. One of the interests to combine these two techniques PAT/DOT is that in time domain they have in common the excitation equipment and will permit to built a compact device. PAT has a high tissue specific optical contrast from light absorption but it does not provide information about the scattering properties of the medium. The fusion of PAT/DOT permits to quantify PAT images thanks to the estimation of the light fluence distribution with DOT [START_REF] Bauer | Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography[END_REF]. DOT can bring information on the scattering coefficient distribution. An integrated PAT/DOT system can take advantages of both DOT providing scattering images and PAT generating high resolution absorption [START_REF] Li | Integrated diffuse optical tomography and photoacoustic tomography: phantom validations[END_REF].

With PA, it is possible to acquire US image if the US transducer has both emission and reception modes then three imaging modalities can be combined.

Other combinations. Recently, other possible combinations with DOT have been considered.

Merging with PET, which is a functional and molecular imaging producing images by recording gamma rays that are emitted from positrons released from radioactively labelled substances, was pursued [START_REF] Liang | Biomedical Optical Imaging Technologies[END_REF]. Functional DOT can be combined with monitoring devices such as electrocardiography (ECG) on an ambulatory system for muscle physiological monitoring during exercise applications [START_REF] Hu | Ambulatory diffuse optical tomography and multimodality physiological monitoring system for muscle and exercise applications[END_REF].

Research on multimodal imaging are more and more exploited and for optical imaging is a big step to go to clinical. A commercial system developed by Bruker in vivo Xtreme combines multi-imaging technologies but it is for small animal imaging.

To conclude, hybrid imaging has great advantages for DOT to overcome its resolution limitation by using the spatial information derived from other imaging modalities to guide the inverse optical problem and to provide complementary information to improve the diagnosis. However, multimodality imaging has also drawbacks which reduced some potentials of DOT when it is used on its own: the cost which is increased by the use of several modalities, the volume, potentially the ionizing, the technological constrains for the compatibility between the combining modalities which limit some freedom instrumentation design and the complexity to merge all the data into the data processing. problem explanations are developed with more details than in the previous chapter. In a second section, the various multispectral methods and modifications of the algorithm proposed are explained. In a third section, we describe the data analysis to assess the detection and quantitation performances. In a fourth section, we evaluate with numerical investigations the robustness and reliability of our algorithm to reconstruct, to detect and to quantify scattering tissues. A first simulation study examines the noise stability and compares the different algorithmic modifications. This experiment permits to select the optimal algorithm for the rest of the manuscript. A second simulation study evaluates the various multispectral methods where the optimal one is retained for the Chapter 4.

Direct and inverse problems

In order to retrieve the optical variation of a medium in a tomographic representation, we use a reconstruction method implying a direct model (or forward model) and an inverse problem. In this section, we will explain in details each block of the reconstruction process.

Forward model

Diffusion approximation

Due to the high complexity of biological tissues, they are very diffusive. In this condition, the diffusion approximation (DA) is a valid model of propagation as long as the acquisition geometry respects the distance condition (higher than 1/µ ′ s ) and is commonly used as explained in Chapter 1. In our reconstruction method, we employ the DA as forward model. At each space position r and time t, the diffusion equation links the photon density ϕ( r, t) to the source S( r, t) and to the optical properties of the traveling medium through the diffusion and absorption spatial coefficients distributions D( r) and µ a ( r) as formulated in Eq. 1.5. If S( r, t) is a Dirac source δ(t) * δ( r -r s ) at position r s and time t = 0, then the photon density ϕ( r, t) is called the Green's function G( r s , r, t). The differential equation 1.5 becomes:

1 c ∂G( r s , r, t) ∂t -∇.(D( r) ∇G( r s , r, t)) + µ a ( r)G( r s , r, t) = δ(t) * δ( r -r s ) (2.1)
For the rest of the mathematical development, we will work on this equation 2.1.

Note that if the medium is illuminated by a collimated source, it is equivalent to an isotropic source S( r, t) at 1/(µ ′ s ) within the medium. The model must be supplemented with boundary conditions that express the reflections of light at the interfaces of the medium; it implements that there is no light coming from outside (expect from the source illumination which is handled by S). To take into account these constraints, we apply the modified Robin boundary condition [3,[START_REF] Martelli | Light propagation through biological tissue and other diffusive media theory, solutions, and software[END_REF][START_REF] Montcel | Tomographie optique diffuse résolue en temps: Applications fonctionnelles en neurosciences[END_REF][START_REF] Schweiger | The finite element method for the propagation of light in scattering media: Boundary and source conditions[END_REF]:

αD ∇ϕ( r, t). n out = ϕ( r, t) (2.2)
with α = 2 R+1 R-1 is a coefficient taking into account the index mismatch between the medium and the outside where R is a parameter governing the internal reflection at the boundary (being the Fresnel reflection coefficient between two media) and n out is the outwards-oriented normal vector.

The system of equations 2.1 or 2.2 is to define photon densities ϕ or G of the differential diffusion equation when µ a , D and the geometry are known.

DA solutions

Different methods can be used to solve the DA depending on the complexity and geometry of the medium. For some simple geometries (infinite, planar, cylindrical or slabs), analytical solutions of DA exist [START_REF] Arridge | The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis[END_REF][START_REF] Jacques | Tutorial on diffuse light transport[END_REF]. For example, in the case of an infinite homogeneous medium, we have this analytical formula:

G( r s , r, t) = c (4πDct) 3/2 e -r 2 4Dct -µact when r s = 0 (2.3)
To solve the diffusion equation in the case of media with arbitrary shapes and/or arbitrary optical parameter distributions, it is necessary to use numerical methods such as finite-element method (FEM) [START_REF] Arridge | Optical imaging in medicine: II. Modelling and reconstruction[END_REF][START_REF] Gao | Optical tomographic mapping of cerebral haemodynamics by means of time-domain detection: methodology and phantom validation[END_REF][START_REF] Schweiger | The finite element method for the propagation of light in scattering media: Boundary and source conditions[END_REF] or finite-volume method (FVM) [START_REF] Hervé | Multi-wavelength and time-domain diffuse optical tomography data processing by using a material basis and Mellin-Laplace transform[END_REF]. These methods permit to spatially discretize the equation by using a finite partitioning in tetrahedra of the medium. We will use FVM, since for regular mesh (as in our case) it is more accurate. The surrounding medium has a µ ′ s of 12 cm -1 and a µ a of 0.07 cm -1 . The presence of the inclusion modifies the Green's function through µ a ( r) and µ ′ s ( r) and the difference with the reference Green's function (homogeneous medium without inclusion) starts at short time (i.e. early photons are absorbed and thus not collected) when the inclusion is at 10 mm depth. At deeper layers such as 20 mm depth, the change starts later because it is later photons that are involved through some of them being absorbed. To better see the effect of the inclusion, the contrasts, defined as the ratio between the difference of the reference and perturbed Green's function are plotted in Figures 2.2 (c) and (d). In both cases, at a given time the contrast increases logarithmically in accordance with the increase of the absorption variation of the inclusions in both case (black arrows). However, the slope of the contrast for inclusions at 20 mm depth raises slower and later than at 10 mm depth.

Multispectral time-resolved measurements Time-resolved measurement model

In order to get 3D maps of material composition, the probed medium is scanned at several TR measurement M sd (t) (in our case a count of photons) for source s and detector d at time t is the result of 3 steps: the light injection in the medium, the light propagation and the detection.

From the mathematical point of view, a TR measurement is the time convolution (we suppose that the source and detector are seen as points):

M sd (t) = S s (t) * G sd (t) * D d (t) (2.4)
where the operator * is the time convolution product. S s (t) is the input signal of the source and D d (t) is the detector response. From now to simplify the notations, the position r s and r d of source and detector will be replaced by the indices s and d, respectively. So G s ( r, t)

is equal to G( r s , r, t), G d ( r, t) is equal to G( r, r d , t) and G sd (t) is equal to G( r s , r d , t). The
Green's function is symmetric; G( r, r ′ , t) is equivalent to G( r ′ , r, t). Equation 2.4 enables us to introduce the instrumental response function (IRF) of the experimental setup which is defined as the combination between the time functions of the source and the detector (in practice not separated). So, we can write:

M sd (t) = IRF sd (t) * G sd (t) (2.5)
The IRF corresponds to some non idealities introduced into the measurements (source and detector delays, widths, tails, decays ...). Figure 2.3 illustrates this equation. In some studies [START_REF] Valim | The effect of temporal impulse response on experimental reduction of photon scatter in time-resolved diffuse optical tomography[END_REF] and in our lab [START_REF] Laidevant | Méthodes optiques résolues en temps pour la tomographie de fluorescence dans les milieux diffusants[END_REF][START_REF] Puszka | Diffuse optical tomography: a time-resolved approach for reflectance measurements at short source-detector separation[END_REF], it was demonstrated that the IRF has an influence on the DOT results. If we want to keep time information from measurements, the significant characteristics of the IRF in the order of importance are (Figure 2.4):

1. Decay: The slope in log-scale of the IRF must be steeper than the slope of the Green's functions which depends mainly on µ a (see Figure 2.2) [START_REF] Laidevant | Méthodes optiques résolues en temps pour la tomographie de fluorescence dans les milieux diffusants[END_REF][START_REF] Torricelli | Time domain functional NIRS imaging for human brain mapping[END_REF] otherwise the information of the medium will be hidden by the instrumental response. The slope of the IRF can be quantified with the exponential decay constant, α defined with the equation of the slope exp(-α × t).

Dynamic range:

The dynamic range is also crucial to obtain information in depth otherwise it will be hidden by the noise at the tail of the IRF.

Full-width at half maximum (FWHM):

The FWHM is the width calculated by the difference between the two extreme values on the curve at which the function reaches 50% of its maximum value. The FWHM of the IRF has to be finer than the FWHM of the Green's function of the studied medium because of the risk to hide information into the Green's function. However currently with the picosecond pulse laser and the current detectors, the problem of the FWHM is mainly solved (see Appendix A).

A full study on the effect of the IRF on the quantification is done in section 3.1 and different experimental IRFs are plotted in Figure 3.2.

The interest and objective of TR signals are to allow to distinguish early and late photons so that the information in both superficial and deep layers can be retrieved. We can see that the early photons transit in the red area which is the shallow layer and the late photons reach deeper layers corresponding to the pink area. Then below the schema, Figure 2.5 shows several experimental temporal distribution of photons which were measured on media with or without a perturbed absorbing object at different depths. These measurements were acquired with FG SPAD. In general, we observe a change of the shape due to the presence of the object and its position in the medium. At 5 mm depth, very early photons are affected by the inclusion and we observe an offset between the reference signal (without the perturbation, red plot) and the measured signal (blue curve) and even the height of the signal is reduced. We have a decrease in amplitude. Then when the inclusion is embedded deeper at 15 mm or 20 mm depth, only the higher order channels are affected i.e. the longer time. And at short time, the two curves are superimposed. At 20 mm depth, it becomes difficult to distinguish the two curves. A perturbation in absorption affects significantly the TR measurements and the depth of the perturbation can be deduced from the time where the measured and reference signals depart.

Multispectral approach

Another dimension is added to the TR measurements for DOT, it is the multispectral approach.

In Chapter 1, we have explained that the main clinical motivation of diffuse bio-optical imaging is that biological tissues are composed of chromophores with specific absorption spectra (Figure 1.1). By probing biological tissues with NIR light, we are able to measure the concentrations and positions of chromophores in the tissue. Multiwavelength (MW) acquisitions in TR DOT will permit to separate contributions of chromophores and determine their absolute concentrations by using the law of mixtures Eq. 1.3 where ε i is chromophore molecular extinction coefficient (L.mol -1 .cm -1 ) which describes the "absorption strength" of the chromophore i. A higher value indicates a stronger absorption.

Mellin-Laplace transform

In order to simplify data processing of time measurements we employ the Mellin-Laplace transform (MLT) [START_REF] Hervé | Time-domain diffuse optical tomography processing by using the Mellin-Laplace transform[END_REF], a temporal integral applicable to the DA equation 2. The general form of the MLT for a time function f is:

M LT (p,n) [f (t)] = f (p,n) = +∞ 0 f (t) (pt) n n! e (-pt) dt (2.6)
By applying the MLT on the diffusion approximation (Eq. 2.1) [START_REF] Hervé | Time-domain diffuse optical tomography processing by using the Mellin-Laplace transform[END_REF], we get:

-∇.(D ∇G (p,n) s ) + (µ a + p c )G (p,n) s = δ( r -r s )δ n + p c G (p,n-1) s (2.7)
The convolution operator ' * ' becomes the discrete convolution product on n:

(f * g) (p,n) = N i=0 f (p,i) × g (p,n-i) (2.8)
The larger the values of p, the more pieces of information are extracted from the diffusion curve per nanosecond. The larger the value of n, the more late "time-windows" are included in the analysis. These two remarks on MLT are illustrated in and p = 3 ns -1 respectively.

Inverse problem

The previous part tells us how to predict measurements for a given distribution of optical properties of the medium, it is the forward model. What we want to do now is to determine 3D maps of material composition from time-resolved optical signals; it is an inverse problem. To do this, differences between actual measurements and predictions are used so as to update the optical properties map. In DOT, it is mainly obtained by using the perturbation theory which gives the variation of the Green's function due to variations of optical properties.

Perturbation approach

The perturbation approach permits to link explicitly the variations of Green's functions G with the variations of optical parameters.

Two common perturbation approaches exist: the Born or the Rytov approximation where a change in µ a and µ ′ s of the heterogeneities is regarded as a perturbation from a reference state [START_REF] Arridge | Optical imaging in medicine: II. Modelling and reconstruction[END_REF][START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF][START_REF] O'leary | Imaging with diffuse photon density waves[END_REF]. For only µ a changing, the Born approximation is the following:

δG sd (t) ≈ - Ω G s ( r, t) * δµ a ( r) * G d ( r, t)d r (2.9)
where δG sd (t) is the small variation of Green's functions between two states of the medium at a time t and δµ a ( r) is the corresponding variation of absorption coefficient. With this Born approximation in TD [START_REF] Fukuzawa | Reduction of image artifacts induced by change in the optode coupling in time-resolved diffuse optical tomography[END_REF][START_REF] Zhang | Combined hemoglobin and fluorescence diffuse optical tomography for breast tumor diagnosis: a pilot study on time-domain methodology[END_REF], a linear system is obtained:

δG sd (t) = Ω W ( r, t)d r δµ a (2.10)
with a matrix W which is called the Jacobian, sensitivity matrix or weight matrix (these different terminologies are used in the literature). The sensitivity matrix is involved in equation 2.9 which represents the effect of absorption perturbations between two states on TR Green's functions on a volume Ω. We only represent the variation in absorption coefficient δµ a with the same fixed reduced scattering coefficient for both states.

The Rytov approximation can also be used to compute variations of G sd . In the literature, this perturbation approach enables to linearize the inverse problem in CW [START_REF] Boas | Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy[END_REF][START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF][START_REF] Kavuri | Sparsity enhanced spatial resolution and depth localization in diffuse optical tomography[END_REF] although with more difficulties in TD [START_REF] Cai | Time-resolved optical diffusion tomographic image reconstruction in highly scattering turbid media[END_REF][START_REF] Namita | CT imaging of diffuse medium by time-resolved measurement of backscattered light[END_REF]. Instead of considering changes in absolute intensity as in the Born approximation, in Rytov approximation changes in logarithmic intensity can be used to derive the sensitivity matrix [START_REF] O'leary | Imaging with diffuse photon density waves[END_REF]. This latter approach is based on this approximation:

δlnf = δf f (2.11)
which is correct for weak variation of an arbitrary function f . Applying on G sd the formula 2.11 and knowing that G 0 sd = G sd -δG sd , it gives the Rytov formulation [START_REF] Arridge | Optical imaging in medicine: II. Modelling and reconstruction[END_REF]:

G 0 sd ln G sd G 0 sd ≈ - Ω G 0 s ( r, t) * δµ a ( r) * G 0 d ( r, t)d r (2.12)
where G 0 sd (t) is the reference Green's function measured at a source-detector couple position at a time t with µ 0 a ( r), and G sd (t) is the Green's function with µ 0 a ( r) + δµ a ( r). So, the Rytov formulation Eq. 2.12 is based on the logarithm ratio of the perturbed signal G sd (t) and the reference signal G 0 sd (t) (development of the Rytov formulation in [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF]).

We tried to implement the Rytov approximation but contrary to the Born approximation we did not succeed in eliminating the IRF from equations so we need to know a priori the IRF to insert it into the code. Our technique to eliminate the IRF with the Born approximation will be explained in section 2.1.2.2. Some tests on a simulation study using Rytov approximation have been done and gave no improvement on quantification and more artefacts in depth appeared on the reconstructed maps (not shown). By contrast Arridge et al [3] and O'Leary et al [START_REF] O'leary | Imaging with diffuse photon density waves[END_REF] have found that Rytov approximation provided improvement and produced more accurate reconstructions on FD measurements. Maybe some future works can be done to have a full comparison study between the Born and Rytov approximations in TR DOT.

Cross iterative Born approximation

Before going to the inverse problem block (Figure 2.1), we will introduce the perturbation approach developed by Dr L. Hervé [START_REF] Hervé | Time-domain diffuse optical tomography processing by using the Mellin-Laplace transform[END_REF] and used in our MW TR reconstruction process. This approach based on the Born approximation called in this report "cross iterative Born approximation" permits to avoid the involvement of the IRF in the processing of the reconstruction algorithm. The developed approach requires a measurement on a reference phantom with known properties, but in practical these properties have to be close to the background of the probed medium. To be in the framework of the perturbation approach, i.e. to have small variations of the absorption coefficient or the concentration, it is recommended to employ an iterative process. Note that we will consider the reduced scattering coefficient µ ′ s known and constant in the homogeneous and perturbed medium.

The process is that we acquire two sets M Aλ sd (t) and M Bλ sd (t) of multi-wavelength (indexed by λ) time-domain measurements on a known reference medium A and on the unknown chromophore composition medium B for a collection of source (indexed by s) and detector (indexed by d)

positions (Illustration Figure 2.5(a)). M Aλ sd (t) = IRF λ (t) * G Aλ sd (t) (2.13) M Bλ sd (t) = IRF λ (t) * G Bλ sd (t) (2.14)
G Bλ sd (t) is unknown in equation 2.14 since µ λ a is unknown.

These two equations are involved in the perturbation approach in the framework of an iterative process. Here we use our cross iterative Born approximation by performing the operation

(2.14) * G Aλ sd (t)-(2.13) * G Bλ(k) sd
(t) and we get at iteration k:

M Bλ sd (t) * G Aλ sd (t) -M Aλ sd (t) * G Bλ(k) sd (t) = M Aλ sd (t) * G Bλ sd (t) -M Aλ sd (t) * G Bλ(k) sd (t) (2.15) = M Aλ sd (t) * δG Bλ sd (t) (2.16) with δG Bλ sd (t) = G Bλ sd (t) -G Bλ(k) sd
(t) corresponding to the variation of the Green's functions between the experiment and the prediction at iteration k. δµ

λ(k) a is equal to µ Bλ a -µ Bλ(k) a . k is the iteration parameter. G Bλ(k) sd (t)
is the estimate Green's function of the medium B at iteration k. According to the Born approximation, we get:

δG Bλ sd (t) = - Ω G Bλ s (t) * δµ λ(k) a * G Bλ(k) d (t)d r (2.17)
which is deduced from the effect of absorption perturbation on TR Green's functions on a volume Ω. So, we obtain the following system:

Y λ(k) sd = -M Aλ sd (t) * Ω G Bλ(k) s (t) * δµ λ(k) a * G Bλ(k) d (t)d r (2.18) Y λ(k) sd = M Bλ sd (t) * G Aλ sd (t) -M Aλ sd (t) * G Bλ(k) sd (t) (2.19)
where Y is defined as the matrix of the combination of the MW TR measurements and/or the calculated measurements with the direct model and it will be used to set our inverse problem.

Implementation of the inverse problem equation

Now that we have a linear inverse problem, we prepare the equation of the inverse problem to be solved thanks to mathematical tools.

Starting from Eq. 2.18, we want to determine the absorption variations. To implement this problem, instead of time discretizing Eq. 2.18 with time steps dt which would be computationally intense, we propose to calculate its MLT (Eq. 2.6). By applying the time convolution properties of the MLT (Eq. 2.8), it gives us the following linear system:

Y λ(p,n)(k) sd = - i+j+l=n M Aλ(p,i) sd,m × m G Bλ(p,j)(k) s,m × δµ λ(k) a,m × G Bλ(p,l)(k) d,m × V m (2.20) Y λ(p,n)(k) sd = m W λ(p,n)(k) sd,m × δµ λ(k) a,m (2.21)
with

Y λ(p,n)(k) sd = e+f =n M Bλ(p,e) sd,m × G Aλ(p,f ) sd,m -M Aλ(p,e) sd,m × G Bλ(p,f )(k) sd,m G Aλ(p,n) sd and G Bλ(p,n)(k) sd
are the MLT of Green's functions of the diffusion equation between the source s and the detector d. V m is the volume of the mesh index m.

Equations 2.20 and 2.21 are linear. They can formally be written as:

Y = W [µ] (2.22)
where Y is the matrix of the combination of the measurements and/or the calculated measurements with the forward model, [µ] is the matrix of the unknowns containing the variations in absorption coefficients δµ a (µ ′ s is supposed constant) and W is the model. W is defined as the sensitivity matrix, Jacobian matrix or weight matrix.

Thanks to the mixture law (Eq. 1.3), our system to solve (Eq. 2.22) can also be written as:

Y = W S[C] (2.23)
where the spectra matrix S is the specific absorption spectra or the molar extinction coefficient spectra of the chromophores with dimensions N λ × N c with N c the number of chromophores indexed by c and N λ the number of wavelengths; then the matrix [C] represents the variation in proportions or concentrations of components in the medium (δC) respectively.

In general form, our system becomes:

Y = W P X (2.24)
where P is a transfer matrix adapted to the nature of X. When we want to work on absorption coefficient variation (δµ a ) called the Indirect method, X is equal to [µ] and P is transformed to a identity matrix (Id). Otherwise when we want to retrieve directly the concentration variation (δC) called the Direct method, X is equal to [C] and P is equal to the spectra matrix S. We will see further the Indirect method and Direct method in section 2.2.1. We will use this equation 2.24 for the rest of the manuscript.

Computation of the update

With our implemented system (Eq. 2.24), we can go to the last step of our reconstruction algorithm: the computation of the update. Firstly, our system Eq. 2.24 is not square and is ill-conditioned. So, there is no exact solution.

We are looking for the closest approximate solution X; the solution is to take the least squares sense. To solve this kind of problem, we use a technique of minimization of a criterion χ 2 with the generalized least squares criterion written as:

χ 2 =|| Γ -1/2 Y (Y -W P X) || 2 (2.25)
A solution is characterized by:

∂χ 2 ∂X = 0 (2.26)
where Γ Y is the covariance matrix of Y . χ 2 represents the difference between the experimental measurements and the model. We note Γ -1/2 Y = L where L will be the left preconditioning matrix which weights the measurements in accordance with their associated noise. L may be equal to Γ

-1/2 Y
in case of generalized least squares. In our case, we choose L as a diagonal matrix with elements 1/σ λ sd in order to whiten the noise where σ λ sd are the standard deviations on Y obtained by deriving the initial Poisson noise of the time-resolved measurements. This choice of L can be seen as the weighted least squares (explained in detail in section 2.2.2).

To relieve the process to solve this ill-conditioned system, we decompose our system by using a transfer matrix R which we will call the right preconditioning matrix. R permits a change of basis. So, the equation 2.25 is transformed to:

χ 2 =|| L(Y -W RX ′ ) || 2 (2.27) with P X = RX ′ (2.28)
And by setting Y ′ = LY and W ′ = LW R, we obtain this following new problem which is an ordinary least squares problem.

χ 2 =|| (Y ′ -W ′ X ′ ) || 2 (2.29)
L and R are called respectively the left and right preconditioning matrices. R may be composed of different matrices which have different roles:

❼ matrix P : the matrix R contain P which can be equal to the spectra matrix S when we want to retrieve the concentration maps. With one wavelength, we want to recover the absorption coefficient and thus P is the identity matrix (Id).

❼ matrix D: R can be composed of the matrix D called the Dictionary matrix. So, we can write X = DX ′ . The two possible functions of this matrix are to weight the voxels with coefficients which may permit to assign different weights on each voxel relative to the distance between the voxel and source-detector probe and/or to weight the voxels spatially by constraining the reconstruction on wished areas. Each column gathers voxels together into areas which are called "words" if desired. With this matrix, reconstruction by including a priori the position and the size of a region of interest is possible and called in this manuscript the spatially constrained method or a priori code detailed in section 2.2.3.3. It is also possible to do reconstruction by layers.

Different forms of R are presented in section 2.2.3 and tested in a simulation study in section 2.4.2.

For example, in the case of P = S and so X = [C], we get this equation:

[µ] λ m = c P λ c X m,c (2.30)
= c P λ c D D m,D X ′ D,c (2.31)
where c the index of chromophore (number of chromophore N c), D indicates the "words" (i.e.

the areas) and m is the mesh index.

In order to determine X, we first have to deal with this new inverse problem Y ′ = W ′ X ′ . We minimize the associated least squares criterion χ 2 (Eq. 2.29) with a method as the conjugate gradient method (as in our case with 5 iterations) to obtain an estimation of X ′ .

Finally, with the matrix multiplication 2.28 the final X can be determined. For large variations between the reference medium and the perturbed medium, we iterate the process as summarized in the diagram in figure 2.1. We inject in the direct model an update of X.

If X = δµ (k) a,m , we update δµ (k) a,m = δµ (k-1)
a,m + γX with γ = 0.5 (typical relaxation parameter). Then we go back to the perturbation approach and the retro-projection with the inverse problem until reaching the desired number of loops or until convergence between the model and the experiment. Most of the time, we fix the stop condition at 10 iterations; it is enough to reach the convergence.

Developed multispectral algorithm

Our developed time-resolved multispectral algorithm will be summarized in this section with two possible multispectral methods. This algorithm is based on the TR DOT algorithm developed by L. Hervé at the laboratory LISA of CEA-LETI presented above. This algorithm can be used for any geometry of the medium and configuration of sources and detectors (for example: cylindrical, spherical, and planar), in 2D and 3D. Some proposal of improvements and adaptations of the code will be given.

Two multispectral methods

As a reminder, our algorithm is based on the diffusion equation Eq. 2.1 which is MLT solved and to linearize the inverse problem we use the "cross iterative Born approximation" perturbation approach.

To determine the concentrations of chromophore with the multispectral approach, we must change the basis of the unknown from δµ a to δC the variation in concentration. We use the law of mixtures:

[µ] = S[C] (2.32)
where S is always the spectra matrix. Depending on the method, this matrix equation is at the beginning or the end of the algorithm. The system to solve is still Eq. 2.24 where P and X can be changed.

With the multispectral approach, we add another dimension to our problem. We are now working with 6 dimension matrices (m mesh, s source and d detector positions, n order, λ wavelength and c chromophore). So, we have to handle carefully the matrices to avoid any errors.

Two algorithms can result from this multispectral algorithm and can be implemented in the reconstruction process: Algorithm 1 called the Indirect algorithm and Algorithm 2 called the Direct algorithm which are described below.

Algorithm 1 Indirect P = Id and X = [µ] loop λ loop iteration k Y ′λ = W ′λ X ′λ ① X λ = DX ′λ δµ λ a ← δµ λ(k-1) a + γX λ ② end loop Final µ λ a end loop C = S -1 µ a ③ Final C Algorithm 2 Direct P = S and X = [C] loop iteration k Y ′ = W ′ X ′ ① X = DX ′ C ← C (k-1) + γX ② end loop Final C 2.2.1.

Indirect algorithm

It works in two steps: first step (① and ②) is the reconstruction of each absorption coefficient at each wavelength on each point of the mesh µ λ a,m independently and second step (③) we switch in concentration thanks to the law of mixtures (Eq. 2.32). So for the first step, we set the P matrix equal to the identity matrix (Id), thus X is equal to [µ] and the inverse problem can be written as ①. At each iteration k, the absorption coefficient is updated and back-projected in the direct model (②). To summarize, the Indirect algorithm is an algorithm wherein the retrieved absorption distributions at each wavelength are first obtained and are then decomposed into chromophore concentrations images thanks to the known spectra of the chromophores S.

Direct algorithm

The Direct algorithm is also called the decomposition in material basis method [START_REF] Zouaoui | In-depth quantification by using multispectral time-resolved diffuse optical tomography[END_REF]. It consists in directly reconstructing the proportions or concentrations of components in the medium without passing through the absorption coefficient. So in equation 2.24 the matrix P is set as the spectra of the absorption coefficients or the molar extinction coefficients of each component c depending on the wavelength λ (matrix S dimension

N λ × N c). We solve directly X ′ containing δC (k)
c,m and thus the inverse problem is defined as ①. This equation is directly implemented in the TR DOT algorithm. Then we update (②) at each iteration k the concentration C with the solution of the inverse problem converted into X which is the differential concentration distribution on the original mesh δC

(k) c,m .
To conclude, the advantage of the Indirect algorithm is that we can study more easily the different combination of the wavelengths (less time consuming) and we can also study the absorption coefficient quantification. We avoid possible problems of ponderation of some wavelengths because for some wavelengths the number of collected photons can be different and so we can have more weight on some data. This depends on the stability of the laser with the wavelength and the response of medium at the light. In the section 2.4.3, the evaluation of these two multispectral methods on a simulation study is developed. In the other chapters, all the results are obtained by using the Indirect algorithm.

Noise consideration

We want to dwell on the methods of how to take into account the noise of the measurements.

The noise consideration is something important to solve correctly the inverse problem although it is often neglected. We already mentioned that we weight the measurements in accordance with their associated noise (L matrix) using weighted least squares. Below we have listed some existing practices to take into account the noise.

Generalized least squares (GLS)

To solve the inverse problem and take into account the noise, the ideal is to apply the generalized least squares (GLS) defined by :

χ 2 =|| Γ -1/2 Y (Y -W X) || 2
(2.33)

χ 2 = t (Y -W X)Γ -1/2 Y (Y -W X) (2.34)
where Γ Y is covariance matrix of de Y .

Simplification of GLS : weighted least squares

We will use a simplification of GLS often called weighted least squares. The definition in wikipedia is "A special case of generalized least squares called weighted least squares occurs when all the off-diagonal entries of Ω (the correlation matrix of the residuals) are null; the variances of the observations (along the covariance matrix diagonal) may still be unequal (heteroskedasticity)."

In our case, this means that we modify the real Γ Y by only taking into account the diagonal elements i.e. only the variances var(y i,i ) and no covariances cov(y i , y j ). So the advantage is that now the inverse of Γ Y is trivial because the modified matrix is diagonal.

χ 2 = i (Y i -(W X) i ) 2 σ 2 Y i (2.35)
We take into account the variances calculated by the Poisson distribution or fluctuations estimated from the variability of the experimental data but considering that the temporal measures (M A and M B ) are independent so there are no covariance elements. This means that Γ M A , the variance matrix of the measurements M A and Γ M B , the variance matrix of the measurements M B are diagonals. Passing through the operations, the MLT of the temporal measurements and Eq. 2.19, which give Y we can deduce Γ Y . And for our weighted least squares approach, we take only the diagonal elements, σ 2 Y i (σ Y i is the standard deviation on Y i with i the diagonal element index).

Other simplifications

Most of the researches on DOT use radical simplications of the GLS by removing the influence of the variance which had to weight the measurements.

It is recommended to avoid using these kinds of simplification because information is not correctly weighted. In this case, it is considering that the variance is equal to 1 for all the measurements and thus Γ Y is an identity matrix which is not true. So it obtains a least squares form, where no noise is taken into account and no variation of the variance.

χ 2 =|| Y -W X || 2 (2.36)
For our algorithm, we will apply the weighted least squares which is a good compromise and integrates well the photonic noise; it is confirmed in the simulation study in section 2.4.2.

Proposed code adaptations

Thanks to our right preconditioning matrix R, we have the possibilities to test several configurations to improve the code. We present below three configurations of R which are tested in section 2.4.2 to evaluate possible improvements in reconstruction and quantification.

Standard code

The first configuration code is the one developed above. We will call this previous multispectral code the standard code where L (left) and R (right) preconditioning matrices are defined as follows: L is still equal to a diagonal matrix with diagonal elements 1/σ

λ(p,n) sd
and R is composed of two matrices: P and D where we add only weighting at each voxels by a coefficient equal to d 2 m , the square of the distance d m between the voxel m and the set of source and detector locations.

Then below we make modifications of the R matrix and keep the same L matrix. These different codes will be tested in section 2.4.2 to evaluate possible improvements in reconstruction and quantification.

R1 code

The second configuration is called R1 code. With this code, we want to keep the sensitivity of the mesh obtained at the first iteration because it is when we get the highest sensitivity in depth. Indeed, we have observed that the sensitivity volume in depth decreases with the iterations. For that, we decide to define a right preconditioning matrix R composed of P and a matrix R1. R1 is a diagonal matrix which we choose in order to keep the sensitivity of the medium at the first iteration J (0) m at all the other iterations. J the sensitivity volume at iteration k is defined as:

J (k) m = sdn W (p,n)(k) sd,m (2.37) We set R1 (k) m = J (0) m J (k) m .
So, it involves weighting the voxels m with the sensitivity volume of the first iteration. This algorithm follows the same process as the FOCUSS algorithm [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: A re-weighted minimum norm algorithm[END_REF]. Additionally, we weight the voxels with elements equal to d m , the distance between the voxel m and the source-detector probe location or d 2 m . We define the code included d m in its weighting elements, the R1 1 code. And otherwise with d 2 m , we call it the R1 2 code.

Spatially constrained method

The third configuration is called a priori code or spatially constrained method. By incorporating geometric constraints in the reconstruction, it was demonstrated in a review [START_REF] Pogue | Implicit and explicit prior information in near-infrared spectral imaging: accuracy, quantification and diagnostic value[END_REF] that prior information in NIRS improves the accuracy for retrieving the expected optical parameters and also in FD DOT [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF]. The Dictionary matrix D can be modified to add a priori geometry by weighting only the voxels of the a priori object. For example, with a system Y = W X to solve and X with the form

X = x ′ 1 .D 1 + x ′ 2 .D 2 (x ′ 1 and
x ′ 2 are unknowns), we can solve the system on X ′ :

Y = W [D 1 D 2 ]X ′ (2.38) with X = [D 1 D 2 ]X ′ (2.39)
We constrain the reconstruction on an expected area at the correct depth where the probed object should be. We use two "words" in the Dictionary matrix and we make 2 conjugate gradient iterations. For this a priori code, the D matrix is composed of two columns (D 1 and D 2 corresponding to zone ① and zone ②) : the first column is composed of all its elements equal to 1 for the background (zone ①) and the second column for the object (zone ②) is composed of 1 when the voxels are inside the expected area and 0 elsewhere. With the same approach, it is also possible to reconstruct layer by layer.

Because of the low resolution of DOT, one of the issues is to do multimodal imaging and to use a high resolution imaging like MRI to bring spatial information seen as a priori for DOT.

Moreover the interest of constraining the reconstruction in FD DOT to improve the quantification was confirmed in Ref. [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF]. That's why we are also exploring a priori approach in TD and see if the quantification can be improved with this kind of algorithm.

Future possible improvements. Some researches have been done on sparsity reconstruction because image reconstruction of DOT appears in a sparse format. The idea is to incorporate this sparseness prior into the image regularization formulation to improve the reconstructed resolution [START_REF] Cao | Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm[END_REF][START_REF] Okawa | Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization[END_REF]. This approach is not developed in this manuscript.

Quantification data analysis

In order to assess the performances of the results obtained from the developed algorithms, a protocol of parameters analysis has to be defined.

For diffuse optics instruments, standardized protocols exist thanks to some successful attempts for performance assessment which have been pursued in the last ten years by reaching consensus in large inter-laboratory studies. As for example, the MEDPHOT protocol [START_REF] Pifferi | Performance assessment of photon migration instruments: the MEDPHOT protocol[END_REF] was conducted for diffuse spectroscopy of homogeneous media, the BIP protocol [START_REF] Wabnitz | Performance assessment of time-domain optical brain imagers, part 1: basic instrumental performance protocol[END_REF] was set for basic performance assessment of time-resolved systems and the nEUROPt protocol [START_REF] Wabnitz | Performance assessment of time-domain optical brain imagers, part 2: nEUROPt protocol[END_REF][START_REF] Wabnitz | Performance assessment of time-domain optical brain imagers: a multilaboratory study[END_REF] was done for dealing with imaging in heterogeneous diffusive media. Since the nEUROPt protocol addressed in general imaging systems, some of its assessment parameters are employed in this manuscript and additional tests are developed. More specifically, the nEUROPt protocol assesses 3 key features that are sensitivity to localized small absorption changes (contrast and contrast-tonoise ratio), spatial resolution (depth selectivity and lateral resolution) and quantification of absorption changes (accuracy and linearity).

In the following, different parameters for signal and data analysis are explained and some of them are extracted from the nEUROPt protocol.

Signal analysis

Concerning signal analysis, we can assess the sensitivity of the instrumentation setup to distinguish an object in a scattering medium. For that, different tools of analysis exist:

Contrast (Const). The contrast is the first parameter which evaluates the effect of variation in absorption between a reference and a measured temporal signals for a given source-detector couple at the same position and at a specific wavelength; M Aλ sd (t) and M Bλ sd (t) respectively.

Const λ sd (t) = M Bλ sd (t) -M Aλ sd (t) M Aλ sd (t) (2.

40)

Signal-to-noise ratio (SNR). The signal-to-noise ratio (SNR) compares the level of a desired signal to the level of noise.

SN R

Aλ sd (t) = M Aλ sd (t) σ Aλ sd (t) (2.41)
where σ Aλ sd (t) is the standard deviation defined as the square of the variance of the TR signal M Aλ sd (t). The highest is SNR, the better information are "visible" and can be extracted. SNR can be increasing at late photons by using gated detectors.

Contrast-to-noise ratio (CNR). The contrast-to-noise ratio (CNR) affects the detection of a small variation in absorption coefficient. The CNR is the multiplication of the contrast and the SNR.

CN R λ sd (t) = Const λ sd (t) × SN R Aλ sd (t) (2.42)
CNR can be increasing at late time by using gated detectors. Indeed, exploiting the gated acquisition enables to collect photons in a selected temporal range for instance after the peak of early photons and thus it increases the overall number of counts recorded in that temporal region which improves the CNR.

Data analysis

The data analysis of the reconstruction 3D maps is divided into two classes. First the ability to detect precisely and the capacity to quantify (accuracy and linearity). A visual evaluation of the quality of the reconstruction in DOT is also important and obtained by looking at the appearances and the scale of the reconstructed 3D images.

Detection analysis

About the detection, we assess the localization of the centroid of the inclusion and spatial resolution.

Centroid localization. To determine the centroid of the detected object, we use the Matlab function regionprops in Image Toolbox on a 3D defined support. So we get the coordinates x, y and z (depth) of the inclusion. In our case, the support volume was determined as the values above the threshold of 50% of the difference between the maximal absorption or concentration perturbation and the background.

Spatial resolution. The spatial resolution is quantified using the FWHM of the z-or y-profiles crossing the center of the inclusion extracted from the 3D maps.

Quantification analysis

The quantification of the composition of biological tissues is a crucial feature for diagnosis. Accuracy and linearity are important because they are directly related to the ability of quantifying chromophore concentrations and their variations.

Accuracy. We evaluate the accuracy on the retrieval of ∆µ a or ∆C with a relative error ǫ defined as (for absorption variation)

ǫ i = ∆µ recovered a,i -∆µ expected a,i ∆µ expected a,i (2.43)
Different ways exist to calculate ∆µ recovered a,i of each region i. It can be the maximal variation in a studied region, written ∆µ max a,i which is the difference between the maximal absolute value and the background value. In clinical circumstance, the quantified results expressing the maximal concentrations are more evident and comprehensible parameters for clinicians. However, because of the low resolution of DOT, one of the solutions for obtaining the retrieved reconstructed value is to calculate an average of the reconstructed ∆µ vol a,i or ∆Ci vol by taking the integral over a certain volume and by subtracting the background value. The background absorption coefficient or concentration is determined by taking the mean in an area without perturbation of the inclusion. The wished volume can be defined by a threshold of a percentage of the maximal value or a spatial region. In this manuscript, by knowing a priori the theoretical position and shape of the inclusion, we often compute the average over the a priori volume and compare it with the maximal recovered value.

Linearity. The term "linearity" is seen in this manuscript as an affine function of the form ax + b. A simple linear regression can be used to assess linearity. Linear regression consists in finding the best-fitting straight line through the points. A common way to summarize how well a linear regression model fits the data is via the coefficient of determination, R 2 . This can be calculated as the square of the correlation between the observed values and the predicted values. If the predictions are close to the actual values, we would expect R 2 to be close to 1.

The coefficient of determination ranges from 0 to 1. We characterized the linearity by the slope and the degree of freedom adjusted R 2 of the linear fit.

For more precision in section 3.2.2.5, we evaluate the linearity on the retrieval of ∆µ a by fitting the dependence of ∆µ recovered a,i vs. ∆µ expected a,i using a 2 nd order polynomial and looking at the linear and non-linear terms.

Evaluation of the MW TR DOT reconstruction algorithm with simulation studies

In this section, we evaluate with numerical investigations the robustness and reliability of our system to reconstruct, detect precisely and quantify objects in scattering tissues and test the various possibilities of the reconstruction program. Two simulation studies were investigated: the first one is to verify the robustness of the algorithm and estimates the best adaptation codes for detection and quantification in depth. The second one is to evaluate the optimal multispectral algorithms on a gradual series of inclusions composed of various concentrations of Hb and HbO 2 . But first we will explain the simulation process.

Simulation process

To assess the algorithms that will be used together with experimental investigations shown in the next chapters, a simulation tool has been developed based on simulated measurements impaired by non-idealities of IRF and statistical noise representative of experiments. This tool permits to better understand the physical mechanisms of TR DOT concerning the quantification in depth and to rapidly test various geometry configurations and various phantom compositions.

The simulation procedure involves first to define the desired simulated phantom composition, source-detector probe and scan geometries. Then to generate the Green's functions through this configuration, we solve the equation diffusion with our MLT approach. As explained in [START_REF] Hervé | Time-domain diffuse optical tomography processing by using the Mellin-Laplace transform[END_REF], by increasing the value of p of the MLT the value of the TPSF is approaching the analytical solution.

So we use for the simulation tool, p = 1/dt where dt = 0.0041 ns and n = 1. The next step is the time convolution of the theoretical or experimental IRF with the Green's functions. Then each curve is multiplied by a scaling factor to get the desired photon integral and the dark count rate (DCR) is added. DCR is the number of counts (uncorrelated with laser pulses) recorded for each time bin (channel) of the TCPSC board depending on the IRF. The DCR affects the dynamic range of the measurements and thus limits the detection in depth. The DCR of the detector is associated with background fluctuation. Finally, photonic noise described by Poisson noise (the main noise contribution) is added to measurements. Figure 2.7 summarizes these steps. Noisy simulation measurements can then be processed via the reconstruction algorithm the same way as real experimental measurements detailed before and give 2D or 3D reconstruction results. To obtain the best performance achievable with our TR DOT reconstruction system, for simulation studies, we apply an "ideal" IRF (close to a Dirac impulsion); in our case it has a FWHM of 4.1 ps (time step of the TR measurements) and an amplitude of 1. We also want to test different noise conditions to have statistical noise analysis. We decided to consider Poisson noise and for the second simulation investigation we include also the background fluctuation based on the DCR. Studies were carried out with a total number of counted photons of 10 7 for each measurement.

Comparison of the different proposed codes on one-wavelength numerical experiments

Objectives and framework

In this simulation study, we want to evaluate different aspects of the reconstruction algorithm.

The main purpose is the determination of the best adaptations of the code. The other questions are the effect of Poisson noise, the detection and quantification capacity, the robustness and viability of the reconstruction process and the effect of source-detector distance.

To answer these concerns, we decided to work on a simulated phantom (see Figure 2.8 (a)) composed of a homogeneous background with µ a = 0.07 cm -1 and µ ′ s = 12 cm -1 without or with a 1 cm 3 spherical inclusion of ∆µ a = 0.15 cm -1 embedded at 3 different depths (10-15-20 mm). The depth is defined as the distance between the surface and the centroid of the inclusion. We simulate a "large" scan with 30 source positions and for each source, two detectors at 15 or 30 mm source-detector distance (ρ) (see Figure 2.8 (a)). The full volume was 80×66×36 Id each voxel

d 2 m standard 1 Id each voxel d m R1 2 code Id each voxel d 2 m R1 1 code Id each voxel d m a priori code Id ① background - ② a priori inclusion Table 2
.1: Summary of the 5 proposed codes evaluate in this simulation study. All the terms are explained in section 2.2.3. R is the preconditioning matrix composed of P which can be the spectra matrix S in multispectral approach and D the dictionary matrix. D has the role to change the spatial basis of the reconstruction with "words" and/or to weight the voxels with a coefficient. d m the distance between the voxel m and the source-detector positions.

To study the effect of Poisson noise and the sensitivity, we did 10 noise draws in the simulation process to generate simulated measurements for 10 numerical experiments. To avoid any other non idealities, we apply a dirac IRF.

Concerning the codes, we evaluate 5 different proposed codes summarized in Table 2.1 and related to explanations in section 2.2.3. The a priori code consists in constraining the inclusion reconstruction on a 1 cm 3 sphere at the expected depth.

This study is restricted to monochromatic measurements on one type of inclusion where the matrix P is equal to the identity matrix (Id).

Simulation results

Figure 2.9: 3D reconstructions of simulation study 1 using standard 2 code. On the top, reconstruction without inclusion and on the bottom with inclusion ∆µ a = 0.15 cm -1 at 10 mm depth.

For the case z = 10 mm, Reconstruction results are presented in 2D slices on Figure 2.9 and 2.10 by employing the standard 2 code. On the profiles shown on Figure 2.11, we compare all the configurations.

Without inclusion. On the 3D reconstructions of Figure 2.9 without inclusion, the noise slightly affects the retrieved µ a distribution by less than 0.003 cm -1 . The area of fluctuation corresponds to the scanned area; indeed a larger inhomogeneous area is obtained with ρ = 30 mm. It seems that at a short interfiber distance, the noise affects more the reconstructed maps (see second row Figure 2.11) because the sensitivity is higher.

With inclusion at different depths. Now with the inclusion, the probed object appears as a spot surrounded by a quite homogeneous background whose µ a is close to the expected value of 0.07 cm -1 . Regardless the volume, the position of the inclusion in both slices (y-z and y-x) is accurately reconstructed (for values see Table 2.2). In addition we clearly see that the retrieved absorption perturbation µ a decreases significantly with the inclusion depth. The expected absolute absorption perturbation equal to 0.22 cm -1 is only obtained for the inclusion at 10 mm depth. Artefacts represented by white spots below and around the absorption spot are probably due to the shape and the coverage of the probe scan which affect the reconstruction in the x-y plane. In the z plan only when the inclusion is at 10 mm depth, these artefacts appear causing by the proximity of the inclusion to the surface. In fact, the configuration of the probe is in L configuration with one source and two detectors one above and the other on the right of the source. So, we can see artefact patterns forming a cross when ρ is equal to 30 mm.

Comparison of each code. The profiles on Figure 2.11 enable to compare the different codes.

The y-profiles are plotted by fixing x = 0 and z at the expected depth. For the z-profiles, we set x = 0 and y = 0. As expected, the profiles of the standard 2 code (red plots) are more square (see y-profiles) and almost reach the expected profiles. For the standard 1 code , the z localization is correct for 10 and 15 mm but it is significantly underestimated for the inclusion at 20 mm. Concerning the quantification of µ a , this latter code gives the less accurate values which are strongly underestimated. We get an overestimation of the inclusion depth with R1 2 code (purple curves) (z values in Table 2.2). Moreover, the profile at 10 mm depth looks like a peak with a lower FWHM than expected and a hollow of the plots between z = 20-30 mm.

With R1 1 code, the shapes of the plots are not respected and the localization at 20 mm depth is wrong. The profiles of the a priori code (gray plots) retrieve the expected profiles for all the depths and thus give an accurate reconstructed µ a . So, the a priori code seems to be the ideal code but the constraint is that spatial a priori information is needed. By analyzing the effect of the photonic noise, the general observation is that it doesn't affect the results. So, in the reconstruction algorithm, we are well taking into account the photonic noise with a weighted least squares approach.

Discussion and conclusion

This first simulation study confirms the robustness and reliability of our system to reconstruct and detect precisely in lateral (x and y) directions taking into account the Poisson noise by weighting it with our L preconditioning matrix in a weighted least squares approach. Simulated a medium without and with an absorbing object affirms the reliability of our reconstruction system to detect small object of 1 cm 3 with a radius of 0.62 cm. The probe configuration provides effects only with close object to the surface where we have a higher sensitivity.

The conclusion about the different codes is that the best code in this kind of configuration is the a priori code, if we have access to spatial information to constrain the reconstruction. This code improves significantly the quantification at all the depths. When the localization is unknown, the standard 2 code is an interesting code for TR DOT but underestimations of µ a are observed for high depths. Weighting the voxels with the square of the distance to the source-detector position improves detection in depth. The standard 2 code allows better shape than R1 2 code in case of reconstruction voxel by voxel. By obtaining a sparse system, the R1 code can probably play a role in the quantification improvement.

To conclude, for the rest of the manuscript we apply the two best codes: the standard 2 code and a priori code. In this study, we only use one inclusion at one wavelength, in the next section results of multicomponents inclusions will be presented.

Multispectral simulation study

Objectives and methods

A multispectral study in TR DOT has been done in [START_REF] Zouaoui | In-depth quantification by using multispectral time-resolved diffuse optical tomography[END_REF] with simulations including an experimental IRF of a gated and non gated SPAD. Only one phantom experiment was tested. To better understand the quantification behavior and the different multispectral algorithms, and to avoid perturbation due to IRF (a study on the effect of IRFs is presented in section 3.1), we decide to do the same simulation study as [START_REF] Zouaoui | In-depth quantification by using multispectral time-resolved diffuse optical tomography[END_REF] but with a "Dirac" impulse of amplitude 1 as IRF.

In this second simulation study, we simulate the scanning of multichromophore (Hb and HbO 2 concentrations) phantoms at three wavelengths 750, 800 and 850 nm. The phantom composition is: a homogeneous background (µ ′ s (λ 750 nm ) = 12 cm -1 constant) with µ a (λ 750 nm ) = 0.12 cm -1 2.12(c). The shape of the inclusion was a cylinder of 8 mm diameter and 22 mm length (see Figure 2.12(a)). The phantom scan was performed by moving the probe composed of one source and two detector positions (ρ = 15 mm) on 6 × 5 positions (distance between two points: 7.5 mm for both directions, see 2.12(a)).

We simulated Green's functions by using the finite volume method on the discretized medium (10478 nodes, 60×50×25 mm 3 medium volume) at three wavelengths (750-800-850 nm) in order to decompose in material the studied medium. The simulated measurements are completed by adding background noise (calculated with DCR) and the photonic noise independently to each channel. To simulate the contribution of the DCR similar of a SPAD detector (DCR = 297 000 counts), we add a offset of 244 photons per time channel to the TPSF. By taking into account these noise sources, the TPSF has a dynamic range of around 4 decades. HbO 2 ((a)), in molar quantities we calculated the total amount of each chromophore inside the inclusion by integrating over the theoretical volume (a priori volume), by taking the integral over an integration support defined with a threshold of 30% of ∆C max (the difference between the maximum reconstructed concentration value and the background reconstructed concentration value) or by integrating over the actual volume of the inclusion with an extension of 3 mm. for 750 nm. Because of the photonic noise and the background noise, the contrast is very noisy for late times and precludes the detection of an object at depths. observed for these two techniques so they seems to be equivalence in simulation.

Multispectral simulation results

In clinical circumstance, the quantified results expressing the maximal concentrations are more evident and comprehensible parameters for clinicians. In case of spatial information with another imaging modality, it allows to measure the integral over the given volume which is equivalent to our quantification technique over the a priori volume of the heterogeneity.

Discussion and conclusion

To evaluate the quantification in depth in DOT and the two different MW methods defined in section 2.2.1, this second simulation study investigated a gradual range of multichromophore phantoms with an inclusion at 10 mm depth. The values chosen for the series of phantoms are in the range of medical interest. Indeed, in literature it was found in invasive breast cancer lesions that HbT is two times higher in malignant lesions (mean maximum 102 µmol/L) than in benign lesions (mean maximum 55 µmol/L) [START_REF] Zhu | Early-Stage Invasive Breast Cancers: Potential Role of Optical Tomography with US Localization in Assisting Diagnosis[END_REF]. This multispectral study demonstrates the ability to separate and quantify chromophore concentrations in depth with our TR MW DOT system from external measurements. Correct localization and good trend of quantification have been obtained. The non-linearity of the quantification of the maximal concentration for high concentration shows the current limit of our MW TR DOT system which will be confirmed with the experiments presented in Chapter 4. Also in [START_REF] Zouaoui | In-depth quantification by using multispectral time-resolved diffuse optical tomography[END_REF], we did a simulation study with the object at 16 mm depth and an important decrease of the retrieved concentrations into the inclusion from 10 to 16 mm depth was observed.

Following the comparison results, the Indirect reconstruction algorithm (Method 1) is retained and will be used in chapter 4. It gives slightly better quantification and it is more flexible about the choice of the wavelengths for the final reconstruction results in chromophore through a material basis.

Some improvements must still be done in the reconstruction algorithm for the multispectral approach to improve the quantification ability on high absorbers and deep objects.

Chapter 3

Detection and quantification of deep absorption perturbations using TR

DOT

The objective of this chapter is to demonstrate through a full study the performances of a whole TR DOT system -the instrumentation and the reconstruction algorithm -to detect and quantify the absorption coefficient distribution of a highly diffusive medium in real conditions.

The issue of quantification in a TR DOT realization in reflectance geometry is addressed.

In this part of the manuscript, we restrict the investigation on the scanning at one wavelength (820 nm) of media composed of one small absorbing object immersed in a scattering medium.

As absorbing inhomogeneities within the homogeneous medium, we use five totally absorbing inclusions (black cylinders) of different volumes. Indeed, in [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 2: experimental implementation[END_REF][START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF] it was demonstrated that a variety of black inclusions with different sizes are equivalent to given absorption perturbations (∆µ a ) over a spherical volume of 1 cm 3 .

This chapter is divided into 3 sections. Before going to the experimental campaigns, we have first exploited the influence of the shape of the instrumental response function (IRF) by examining separately each IRF's characteristics and different medium compositions. We used theoretical IRFs for this study but also real experimental IRFs. We pointed out that the main parameter perturbing the optical information of the medium is the down slope which can hide this information and distort the quantification interpretation if the slope is not steep enough.

The dynamic range is also important to probe in depth. Then two experimental campaigns which were carried out in June 2015 in Milan and additional numerical experiments are presented. The second section investigates the first campaign whose aim was to characterize the potential of TR DOT using single-photon avalanche diodes (SPADs) in non-gated or fast-gated modes (detectors developed by the Politecnico di Milano) to quantify the absorption coefficient in depth at one wavelength. The potential of fast-gated SPAD for detecting in depth up to 25-30 mm has already been demonstrated in [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF]. To better understand the physical mechanisms of TR DOT concerning the quantification in depth, we mimicked all the experimental measurements by numerical measurements simulating equivalent phantoms. We also evaluated the influence of applying the spatially constrained method on real measurements. The second campaign has the same protocol as the previous one but its purpose is the demonstration of the suitability of silicon photomultipliers (SiPMs) for TR DOT. SiPMs have been recently introduced as promising detectors in the field of diffuse optics, in particular thanks to the inherent low cost and large active area [START_REF] Mora | Towards next-generation timedomain diffuse optics for extreme depth penetration and sensitivity[END_REF] and SiPMs represent the most attractive technology for the next-generation of DOT systems [START_REF] Pifferi | New frontiers in time-domain diffuse optics, a review[END_REF]. In this third section, a feasibility study is addressed to evaluate these new promising detectors in the TR acquisition chain to give spatially and quantitative information for TR DOT. This study was also additionally supported by a numerical campaign and a comparison performance study with the SPAD.

Role of the instrumentation on the quantification

The absorption properties of a medium affect the amplitude and mainly the exponential down slope of the Green's functions [START_REF] Laidevant | Méthodes optiques résolues en temps pour la tomographie de fluorescence dans les milieux diffusants[END_REF][START_REF] Torricelli | Time domain functional NIRS imaging for human brain mapping[END_REF]. Thus, to extract information on the absorption distribution in a medium, one of the most important characteristics of the TR signal is the decay (see Figure 2.2). As it was already mentioned in the previous chapter, the TR measurement is the convolution between the IRF of the detector and the Green's function. Consequently, the IRF shape may have an effect on the correct determination of the medium composition.

The shape of the IRF is characterized by a peak defined with three specific parameters (see Figure 2.4) reminded below:

1. The time response quantified by fast up exponential rise and the full width at half maximum (FWHM).

2. The dynamic range associated to the peak and the tail.

3. The decay (also called charge pulse) characterized by a down exponential decay constant, α (ns -1 ) defined in the equation exp(-α × t) describing the decay of the IRF.

Our study focuses on the down slope (called also decay) because of its very likely involvement on accurate quantification and of the lack of previous works dedicated on this issue. The effect of the other parameters has been studied more [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF][START_REF] Valim | The effect of temporal impulse response on experimental reduction of photon scatter in time-resolved diffuse optical tomography[END_REF] but no answer to the problem of quantification has emerged.

Existing detectors for TR DOT and their instrumentation response functions (IRFs)

Before examining of the role of the shape of an IRF, a quick overview of the existing detectors using a time-correlated single photon counting (TCSPC) system for TR DOT and the specific characterizations of their IRFs are detailed here.

First of all, we worked with a TCSPC acquisition chain which enables to acquire at the same time the full TPSF. The principle consists in counting and time tagging the photons one by one.

TCSPC instrumentation converts the signal corresponding to the detected photon into a digital pulse and then measure its time delay relative to a reference pulse. This arrival time is recorded as a count in a histogram, and after multiple laser repetitions builds up into the TPSF (Figure 1.3). The TCSPC board used in the experiment is a SPC-130 (Becker and Hickl, Germany).

The whole procedure for adjusting the settings of the TCSPC board (CFD, TAC and ADC) is provided by the manufacturer in [START_REF] Becker | The bh TCSPC Handbook[END_REF]. We choose a time-to-amplitude converter (TAC) gain of 3 resulting in the smallest time channels possible of 4.1 ps.

Currently, there exists various types of detectors that can be connected to TCSPC board for TR optical measurements listed below.

Photomultipliers (PMs) and Hybrid-photomultipliers (HPMs)

The photomultiplier tube is the first single-photon detector created (in the 1930's). Different technology approaches exist. The most spread option is that an incident photon hits a photocathode which then emits a photoelectron (photoelectric effect). This electron is then accelerated by a high voltage in vacuum and focused on an electron multiplier. This succession of dynodes maintained at increasing potentials enables to multiply the number of electrons (secondary emission) to finally obtain a measurable current at the anode. The involved voltages (photocathode, dynodes) are in the order of magnitude of the kV.

Like for a PM, the hybrid-photomultipliers (HPM) tube involves a photocathode for the photonelectron conversion. The amplification of electrons now takes place in an avalanche photodiode (APD) with a two-step process: first the photoelectron generates a bunch of electrons when hitting the APD and then these secondary electrons are significantly amplified in the avalanche region of the APD. To do this, the photocathode is polarized with a few kV and the APD with a hundred of volts. The principles of PM and HPM are illustrated in [START_REF] Puszka | Diffuse optical tomography: a time-resolved approach for reflectance measurements at short source-detector separation[END_REF].

A comparison study of these two detectors for time-resolved DOT has been carried out in our laboratory in [START_REF] Puszka | Diffuse optical tomography: a time-resolved approach for reflectance measurements at short source-detector separation[END_REF] in 2D. The highest dynamic range and the fast response of HPM enable to obtain contrast close to the ideal case, to obtain better 2D reconstructions and to detect deeper.

The IRFs of both detectors are showed in Figure 3. For the FG SPAD acquisition, we acquire photons using gates opening at different delays in order to record only selected portions of the time-resolved curves (Illustration 3.1(a)) and we adapt the attenuation of the laser emission according to the delay in order to collect the maximum of photons. The first gate was chosen so that it is opened before the first re-emitted photon is collected and closed after the last photons are recorded. In this mode the full distribution of re-emitted photons temporally fits in the gate making this measurement comparable to that performed with a free-running detector (free-running means always active). The non-gated mode (comparable to a free-running detector) and the gated mode are well explained in [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF].

The opening and closing of the SPAD gate were controlled by a delayer synchronized to the laser source and the output of the SPAD module were connected to a TCSPC board. The disadvantage of SPAD is its small active area which affects the collection efficiency of diffused light.

With the gating mode to obtain the full TPSF, some preprocessings have to be done. Our method of reconstitution of the TPSF is to take a window for each delay adapted to a given laser attenuation. Then for each delay, we multiply the signal recorded with the corresponding laser attenuation (Figure 3.1(b)). We sum all the delays together. And to finish we divide by the attenuation of the 1 st delay to recover the true photon counts sent to the medium. Figure 3.1(c) shows the TPSF of non-gated mode where there is no need to build the TPSF and the full TPSF of gated mode reconstituted from different portions. Their corresponding signal-to-noise ratio (SNR) are presented in Figure 3.1(d). For the gated mode, we can see a higher SNR at later time which permits to improve the sensitivity for large time of flight and thus in depth.

Silicon photomultipliers (SiPMs)

SiPM detectors are new detectors with great potential for TR DOT [START_REF] Mora | Towards next-generation timedomain diffuse optics for extreme depth penetration and sensitivity[END_REF]. A SiPM consists of an array of hundreds (or thousands) microcells, each one including a SPAD and its passive quenching resistor. SiPMs feature a time resolution better than 80 ps with a fast tail of just 90 ps decay time-constant. In addition SiPMs bring a number of advantages such as wide active area (up to few square millimeters) that permits to maximize the photon collection, the numerical aperture (both directly reflect into high photon collection efficiency), the high quantum efficiency (easily exceeding 20% at 600 nm wavelength), the broad spectral coverage (350-1100 nm), and the typical compactness of solid-state devices [START_REF] Mora | Fast silicon photomultiplier improves signal harvesting and reduces complexity in time-domain diffuse optics[END_REF][START_REF] Martinenghi | Spectrally Resolved Single-Photon Timing of Silicon Photomultipliers for Time-Domain Diffuse Spectroscopy[END_REF][START_REF] Martinenghi | Time-resolved single-photon detection module based on silicon photomultiplier: A novel building block for time-correlated measurement systems[END_REF]. Essentially, they bring together the advantages of vacuum tubes and microelectronics detectors. SiPMs are also very cheap (much less than 100 ✩ for a single device) because they are used in high quantity for positron emission tomography (PET) imaging. The current limit is the presence of a slow exponential decaying tail (known as "diffusion tail" being due to the diffusion of photogenerated carriers inside the detector) in the IRF of the device (see Figure 3.2(a)). This tail is characterized by a fraction of nanoseconds decaying time-constant starting 2 decades below the response peak [START_REF] Martinenghi | Spectrally Resolved Single-Photon Timing of Silicon Photomultipliers for Time-Domain Diffuse Spectroscopy[END_REF]. This long tail limits the dynamic range of the detector response function to about 2 orders of magnitude, thus possibly affecting DOT reconstructions.

The IRF of each detector presented above were measured in our laboratory or at the Politecnico di Milano where the 5 types of detectors were or are used. These measured experimental IRFs are employed further on in this manuscript (in section 3.1.2.2) for a study on the effect of the IRF on the quantification potential. Their IRF characteristics (FWHM, dynamic range and exponential decay constant, α) are reported in Table 3.1. We remark that SPAD can increase its dynamic range with fast-gated mode.

Investigation of the impact of the IRF

Numerical study with theoretical IRFs

Some studies have investigated the effect on the shape of the IRF on the measurements [START_REF] Valim | The effect of temporal impulse response on experimental reduction of photon scatter in time-resolved diffuse optical tomography[END_REF] and electronical laboratories are still working on the optimization of the IRF shape towards the ideal IRF, a Dirac function. The previous challenge was to detect an object in depth with TR DOT and thus the instrumentation was developed in this direction. For years, the FWHM has been the main parameter improved by the development of new, high speed light sources and detectors. Now, detectors have very fast response which enables the detection of early photons.

Afterwards, the dynamic range was identified as another important parameter to detect in depth so the gating mode was developed by the Politecnico di Milano on the SPAD. Improvements of reconstruction detection with fast gated mode compared to non-gated mode have been demonstrated in [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF]. The technology of photomultipliers has been improved with the HPM which has a higher dynamic range (see Table 3.1). Now, the next challenge is quantification. So, our first step was to evaluate the impact of the characteristics of the IRF on the quantification. Another parameter understudied but also important is the fast exponential decay. As known and reported in Figure 2.2, the absorption coefficient is the dominant parameter affecting the exponential down slope of the Green's function or TPSF [START_REF] Laidevant | Méthodes optiques résolues en temps pour la tomographie de fluorescence dans les milieux diffusants[END_REF][START_REF] Torricelli | Time domain functional NIRS imaging for human brain mapping[END_REF]. Thus the absorption information of a probed medium is principally in this slope. So we decided to study the impact of IRF on this slope by analyzing the contrasts obtained from different absorbing objects.

A simulation study was carried out on various numerical phantoms composed of different absorption perturbations at several depths and different absorption coefficients for the background phantom. The reduced scattering coefficient µ ′ s was fixed at 12 cm -1 . We generated two Green's functions with and without inclusion for each case at a 15 mm source-detector distance. Then, we convolved the simulated Green's functions with the defined IRFs below (the simulation process is explained in section 2.4.1). To study the effect of the down exponential slope of the IRF on the quantification, we simulated multiple theoretical IRF, IRF th , by adding a Dirac pulse at t 1 to an exponential function with an exponential decay constant α (in ns -1 ):

   IRF th (t) = δ(t -t 1 ) where t ∈ [0 t 1 ] IRF th (t) = exp(-α(t -t 1 )) where t ∈ [t 1 + 1 t end ] (3.1)
t is in ns. We arbitrary chose t 1 = 20 × dt with the time step dt = 0.0041 ns. The FWHM is fixed at 0.0041 ns corresponding to dt for all the theoretical IRFs. In order to evaluate one by one parameter, we fixed the same FWHM and we didn't add noise. We tested 7 different numerical IRFs: the Dirac to obtain the optimal contrast and for the 6 others we change the exponential decay constant from 1 to 10 ( α = 1, 2, 3, 4, 5 and 10) (Figure 3.2). The results using these IRFs are indicated by "α k " with k = 1, 2, 3, 4, 5 and 10.

To compare the effect of the theoretical IRFs, we represented in Figures 3.3 and 3.4 the contrasts obtained between the simulated reference (M A sd ) and the perturbed (M B sd ) measurements above the 1 cm 3 spherical absorbing objects defined as:

Contrast(%) = M A sd -M B sd M A sd (3.2)
In Figure 3.3, we used a gradual range of inclusions with absorption perturbations ∆µ a equal to 0.056, 0.087, 0.15, 0.37 and 0.94 cm -1 and the background medium had a fixed absorption coefficient (µ a = 0.07 cm -1 ). In Figure 3.4, we changed the background medium from 0, 0.07, 0.14 and 0.2 cm -1 and fixed the inclusion at 0.15 cm -1 .

Ideal case

In these two Figures 3.

Effect of the IRF decay for various inclusions

The general observation in Figure 3.3 is that the steeper the slope of the IRF (i.e. α high), the higher the contrast and the closer it is compared to the ideal contrast with the Dirac. At the column 10 mm depth, the contrast increases according to the absorption perturbations of the inclusions. For example, for ∆µ a = 0.056 cm -1 and ∆µ a = 0.94 cm -1 , at 5 ns the Dirac contrast is about 32.0% and 87.5% and the α 1 contrast (light blue plots) is equal to 12.0% and 56.6%. This shows that the α 1 contrast becomes closer to the "Dirac" contrast (ideal) for high absorbing objects and thus the effect of the IRF is reduced and the quantification should be more correct. The information contained in the Green's functions are less hidden by the IRF shape. If we now look on the line ∆µ a = 0.15 cm -1 , the depth of the inclusion postpones the start of the increase of the contrast. The information of depth can be detected thanks to this delay. We can observe that a given time like 5 ns, the contrasts decrease with depth and for the "Dirac" contrast the values go from 56.7% to 9.5% from 10 to 30 mm depth. The shape of the IRF is also very important because starting from 20 mm depth, the "α 1 " contrast is nearly equal to 0 and thus the detection is impossible, so as even less the quantification. To detect at 30 mm depth, a longer time acquisition, detectors with a very good sensitivity and IRF with steep with α > 10 are needed.

Effect of the IRF decay for various backgrounds

The change of the background absorption coefficient permits to evaluate how important is the slope of the IRF to detect and quantify heterogeneities in more absorbing environment. On the first row of Figure 3.4 where there is no absorption in the background medium but only scattering, we see that a weak slope of the IRF already reduces the contrast. On the columns of the figure 3.4, we see that this reduction is even more pronounced for high background absorptions. For example, at the column 15 mm depth with µ bkg a equal to 0.2 cm -1 no contrast can be seen for an inclusion with ∆µ a = 0.15 cm -1 when α is smaller than 4 ns -1 . To conclude, the decay of the IRF is an essential parameter to allow accurate detection and quantification and has to be considered. In the reconstruction process, if the IRF is not well taken into account directly or indirectly, it will introduce errors on detection and quantification.

We tested also the influence of the FWHM in Appendix A. In summary, the FWHM introduces only a temporal offset and thus it delays the increase of the contrast.

Effect of the knowing experimental IRFs

After testing the behavior of the contrast with different theoretical IRFs, Figure 3.5 shows the same study as in Figure 3.3 but with 5 experimental IRFs after noise substraction (Figure 3.2):

SiPM, FG SPAD, non gated SPAD, HPM and PM. We decided to subtract the background noise of each detector in order to compare the effect on the shape without its influence. The 4 detectors, SiPM, FG SPAD, non gated SPAD and HPM, have almost the same slope (calculated approximately between 400-800 ps) with α = 7-9 whereas α ≈ 4 for PM. In this latter case, this difference causes contrast to be smaller, less than 10% and for less absorbing objects. Figure 3.5

shows also that at late time the SiPM contrast tends to decrease more rapidly. It is probably due to the "diffusion tail" which hides information at late times. The other difference between each IRF is their dynamic range. In the ascending order, the IRF of SiPM, non gated SPAD, HPM and FGSPAD have around 2.5, 3.5, 5 and 6 decades as dynamic range respectively.

The dynamic range permits to get contrast in depth (i.e. late photons) and thus the possibility to detect deep objects. The dynamic range is also very dependent of the background noise level which can be enhanced with gating mode.

Conclusion and discussion

We demonstrate with this study that the decay of the IRF has an important role on the quantification and on the detection. In the selection of the detector, the exponential decay constant has to be considered in the specifications of a TR DOT setup. In the reconstruction process, a wrong IRF will introduce errors on detection and quantification. Maybe corrections on these errors can be done by knowing a priori the experimental IRF.

In accordance with this study, a nearly ideal detector should have an exponential decay constant superior to 10 ns -1 with a infinite dynamic range and a FWHM equal to 4.1 ps. The current best detectors are FG SPAD and HPM detectors. The SiPM detector follows these two to detect contrast at early time but then it reaches its limits probably due to its "diffusion tail".

Evaluation of the capability of TR DOT to quantify in depth

the absorption coefficient: an experimental and numerical study with single-photon avalanche diodes (SPADs)

Introduction

For this section and section 3. In previous studies, it has already reported results on the sensitivity and spatial resolution of TR DOT using a short source-detector separation scanning scheme [START_REF] Puszka | Spatial resolution in depth for timeresolved diffuse optical tomography using short source-detector separations[END_REF], demonstrating that the adoption of a SPAD has enabled the possibility to detect photons with a long time-of-flight, increasing the sensitivity at depths higher than 20 mm. Further, the use of a fast-gated SPAD allowed to experimentally implement the null-distance approach with advantages in term of contrast, spatial resolution and signal throughput [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF][START_REF] Torricelli | Time-resolved reflectance at null source-detector separation: improving contrast and resolution in diffuse optical imaging[END_REF]. In this section, 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 [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 2: experimental implementation[END_REF][START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF] -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 section is organized as follows: we describe the geometry, the experimental setup, the phantom preparation and the data analysis. We display results using non-gated and gated SPAD for linearity and accuracy both on simulation and phantom experiments. We demonstrate experimentally the benefit of the a priori code on quantification. To finish, we summarize the key findings of the study, address the specific factors that affect TR-DOT quantification, and discuss the implication for specific clinical problems.

Material and methods

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 7.5 mm 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 (Figure 3.6 (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 15 mm (blue crosses) or 30 mm (green crosses) interfiber distance, as represented in Figure 3.6 (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 2 mm (small grey dots in The phantom is in preparation because the tank is half-empty to change the inclusion.

Experimental setup

The experimental setup (whose schematic is reported in Figure 3. m long; Thorlabs Gmbh, Germany). The distance between injection and collection fibers was set to 15 mm and 30 mm depending on the experiment. Photons harvested from each detection fiber were then focused onto a silicon fast-gated SPAD developed by the Politecnico di Milano [START_REF] Boso | Fast-gating of single-photon avalanche diodes with 200 ps transitions and 30 ps timing jitter[END_REF] (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, 3). Schematic experimental setup of the two studies. A light pulse generated with a tens of picosecond pulsed laser at a defined wavelength goes through a variable optical attenuator (VOA) before illuminating a phantom with an absorbing inclusion at a given depth.

Then the backscattered light is collected with two detectors coupled to two TCSPC boards connected to the computer. An electronic synchronization signal at the output of the laser is sent to the TCSPC boards. The green spots correspond to the scan positions of the probe (one source and two detectors).

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 timer 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 10, 15 and 20 mm. 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. [START_REF] Tosi | Fast-gated single-photon counting technique widens dynamic range and speeds up acquisition time in time-resolved measurements[END_REF] 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) . 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 30 mm interfiber distance, we needed 3 delays while for the 15 mm 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 [START_REF] Torricelli | Time-resolved reflectance at null source-detector separation: improving contrast and resolution in diffuse optical imaging[END_REF]. 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 30 mm and 15 mm 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 [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 2: experimental implementation[END_REF][START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF], 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 3.2 shows the dimensions and volumes (EBV) of the objects used in the present study together with their equivalent absorption perturbation (∆µ expected a ) calculated over a 1 cm 3 volume sphere. These equivalent absorption perturbations were determined in [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF] 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 < 10 mm). For this first experimental study with the SPAD, we will work only with the 4 first inclusions i.e. "size 1", "size 2", "size 3" and "size 4". The inclusions were held in a large tank (volume 29×29×14 cm 3 ) through a thin wire attached 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 [START_REF] Spinelli | Determination of reference values for optical properties of liquid phantoms based on Intralipid and India ink[END_REF].

Since multiple measurements are needed to obtain a tomographic reconstruction, during measurements the injection and collection fibers were kept fixed, while the inclusion was moved. This strategy was selected in order to avoid possible artifacts due to the movement of fibers, which were instead bound to the phantom tank using a black PVC plank in which holes for fiber tips were drilled (see Figures 3.6(b) and (c)). Using a motorized 3-axis translation stage, the inclusion was moved over an area of 3.75 × 3 cm 2 at steps of 7.5 mm in both x and y directions (scan area: from -18.7 to 18.8 mm and from -15 mm to 15 mm for x and y direction, respectively). Scans of all the inclusions in the x-y plane were done at different depth (where "depth" refers to the distance from the surface of the phantom to the center of the inclusion) ranging from 10 mm up to 20 mm at step of 5 mm. The depth was automatically changed using the motorized translation stage when a scan on the x-y plane was finished. In total, each inclusion was moved in the three axis x, y and z on 7 × 5 × 4 positions. As reference measurement, for each line, we moved the inclusion far enough (x = -40 mm) to consider negligible its effect.

Preprocessing and reconstruction process Preprocessing

In order to obtain our time-resolved curves used in the reconstruction, some preprocessing are applied on the raw acquisitions. The aim of preprocessing the TPSFs is to remove offsets, select the proper time-window to calculate their Mellin-Laplace transforms (MLT) and reconstruct the TPSFs for the gating mode. We have used the process below for preprocessing TPSFs (Figure 3.1) before calculating their MLT:

1. We took the mean of the repetitions measurements.

2. We calculate the background offset of the signal and subtract it on every signals. We made measurements of the instrumental signal of both boards (board 1 and board 2) without illumination in order to get our background measurements that we subtracted from the experimental signals.

3. We reconstitute the TPSFs for the gating mode as detailed in section 3.1.1.

4. We cut the interesting temporal width of the curve (same for each curve of a scan).

5. We calculate the corresponding variances. Since measurements are affected by photonic noise which is described by a Poisson distribution, the variance of the measurements is equal to the number of counted photons before background subtraction (step 2). For the gating mode, we multiplied the variance of each mean signal for each gate by the square of the associated attenuation coefficient.

For the simulation, we simulated the equivalent absorption spheres of 1 cm 3 and we included the experimental IRF of the SPAD in the simulation process detailed in section 2.4.1.

Reconstruction process

With these preprocessed TPSFs, we can launch the reconstruction process. For the simulation and the experiment, we used the same reconstruction process; the standard code at one-wavelength (see section 2.2.3.1). The stop criterion was set at ten iterations where the convergence is reached.

Spatially constrained method

By incorporating geometric constraints in the reconstruction, we have demonstrated in simulation in section 2.4.2 that spatial prior information maximizes the accuracy in recovery the expected optical parameters. The spatially constrained method (also called a priori code) implemented here consists in using a Dictionary matrix to force the absorbing inclusion on the expected position of the equivalent volume (i.e. a sphere of 1 cm 3 ). We expect that with the a priori code in simulation and experiment the quantification will be improved with such an algorithm.

Data analysis

The objectives of this study are to assess the quantification capacity of TR DOT with SPAD and the effect of depth on this quantification. The analysis of the reconstructed data uses the analysis tools presented in section 2.3 but we apply another method to evaluate the linearity. To be in agreement with the use of the EBV approach, we evaluate the reconstructed ∆µ vol a averaged over a given volume because the maximum value ∆µ max a 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 ∆µ vol a corresponding to a given volume. It is not even simply the product of ∆µ a and volume that is retrieved since the Equivalence Class implies a non-linear relation. Clearly, these equivalence relations are in force only for small objects (e.g. volume 1 cm 3 ) and are ultimately related to the poor spatial resolution of DOT [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF]. In this study, we quantified the absorption variation (∆µ a ) in the DOT reconstruction by comparing the average reconstructed ∆µ vol a,i over a volume of each absorbing object i according to ∆µ expected a,i which is the expected variation. In our case, the average reconstructed ∆µ vol a,i was calculated by taking the average 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 ǫ i defined in section 2.3 with the Eq. 2.43. However for the linearity we apply another method. We evaluate the linearity on the retrieval of ∆µ a by fitting the dependence of ∆µ vol a vs. ∆µ expected a 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:

y = ax 2 + bx + c (3.3)
where y = ∆µ vol a η , x = ∆µ expected a η

, and η = 0.1 cm -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:

slope = 2ax + b (3.4)
Thus, a represents the non-linear distortion that must be referred to the b value. The deviation from linearity can be expressed as:

N L = 2a b ∆x (3.5)
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 session 3.

Results by using gated and non-gated SPADs

For this first experimental study with the SPADs, we will work only with the 4 first inclusions i.e. "size 1", "size 2", "size 3" and "size 4". The results are presented in two parts: first we compare the simulation and experimental results and secondly we see if the quantification can be improved at all the depths with the spatially constrained method.

Comparison between simulation and experiment

The simulation and experimental results obtained with non-gated SPADs of the reconstructed 3D absorption maps are displayed in Figures 3.9 3.2 for equivalent ∆µ a ). The black circle corresponds to a 1 cm 3 sphere and is centered on the expected position of the inclusion.

The absorption perturbation appears as a spot surrounded by a quite homogeneous 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 2 mm is obtained. In Figure 3.9, 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. Figure 3.11 synthesizes the quantified absorption variations ∆µ vol a from the 3D reconstruction absorption results for all the configurations: at 15 or 30 mm interfiber distance, with non-gated or gated SPAD, for simulations and for experiments (dotted and continuous lines, respectively) at 10, 15 or 20 mm depth and for each inclusion. For example, the continuous and dotted brown curves in box "Non Gated" and "ρ = 30 mm" are absorption variations extracted from the 3D reconstruction maps of the absolute absorption distribution of Figure 3.9. 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 ρ = 15 mm. 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 ρ = 15 mm 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 10 mm up to 0.15 cm -1 with relative errors inferior to 30% in the experiment (brown curves in Figure 3.11). Results in simulation for deep inclusion are similar for 15 or 30 mm 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 [START_REF] Quarto | Estimate of tissue composition in malignant and benign breast lesions by time-domain optical mammography[END_REF]. Figure 3.12 displays the percent of relative errors of ∆µ vol a (calculated with Eq. 2.43) between the simulation and the reality and between the experiment and the reality. The same is reported in Figure 3.14 but applying the constrained method. Using the standard code (Figure 3.12), 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 10 mm depth and ǫ = -49% on average with a standard deviation of 4% at 15 mm depth. For inclusion "size 3" at 10 mm, we obtain ǫ = -28% with σ = 6%.

To extract quantitative information on the linearity of the reconstruction, we report in Table 3.3 the coefficients of the polynomial fit of ∆µ vol a vs. ∆µ expected a divided by a reference η = 0.1 cm -1 as defined in section 3.2.2.5. The slope and the N L term are estimated for x = 1 and ∆x = 1, that is around η= 0.1 cm -1 . Regarding 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 = 10, 15, and 20 mm 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 = 10, 15, and 20 mm respectively. This is a strong effect, still not so huge to mask deep changes due to more superficial alterations. The N L coefficient accounts for the non-linearity for increasing ∆µ expected a and is around -20%, at all depths. Thus, the non-linearity with the increased ∆µ expected a is acceptable at least for absorption changes in the order of 0.1 cm -1 . The c coefficient that is the offset at ∆µ expected a = 0 is substantially negligible ( 5%) (Table 3.3).

The coefficients (Table 3.3) 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 ρ = 30 mm is substantially similar to what observed on simulations still with a higher non-linearity (more around 30%). For ρ = 15 mm there is more discrepancy with simulations, clearly due to noisy measurements particularly for the non-gated detector as observed in Figure 3.11. Upon applying the gating, data are slightly more regular, in particular at the larger depth (z = 20 mm).

Results with the spatially constrained method

With the spatially constrained method (Figure 3.13), 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 30 mm interfiber distance, the experimental curves are linear with a small offset. For ρ = 15 mm, 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 ρ = 15 mm than the non-gated mode but no difference is observed for ρ = 30 mm. 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 Figure 3.13.

Thanks to the constrained method (Figure 3.14), 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 Figure 3.14 are almost all positives and for ρ = 30 mm, average and standard deviation (respectively 44% and 14%) of relative errors are lower in absolute and less spreading than in Figure 3.12 (respectively -53% and 22%) without the a priori approach. At the 3 depths, no significant difference is observed between non-gated and gated mode.

Table 3.4 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 ρ = 30 mm. For ρ = 15 mm results are substantially altered with a strong non-linearity and non-systematic alterations. It looks like measurements at ρ = 15 mm 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

We have addressed the quantification of DOT using time-domain reflectance measurements collected with SPAD and using our reconstruction process based on Mellin-Laplace transforms both on simulations and on phantoms. We analyzed the values of the reconstructed absorption perturbation. Five main conclusions can be drawn from the results:

1. The main parameter affecting the correct reconstruction of ∆µ vol a is the depth z. A depthdependence underestimation of the absorption is observed, with a reduction of the slope by a factor of 2, 3 and 6 for z = 10, 15 and 20 mm respectively.

2. The reconstructed ∆µ vol a 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 N L of around 20% for a change in ∆µ a = 0.1 cm -1 (for a 1 cm 3 volume).

3. The adoption of a spatially 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, in [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF] 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 ρ = 15 mm distance are more scrambled and more prone to experimental alterations. The fast-gating does not help here. 5. Experimental measurements on phantoms substantially agree with simulations, at least for the general trends. Apart from some intrinsic higher variability, experimental 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 experimental measurements.

The origin of the depth-and absorption-related underestimation in ∆µ a -already observed in different DOT papers -seems to be indeed due to the general spread of the reconstructed ∆µ a -possibly due to the ill-posedness of the problem -leading to a dilution of the absorption change. Constraining the region of the optical perturbation -as recalled above at item #4completely 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.

New detectors, silicon photomultipliers (SiPMs) for TR

DOT: a feasibility study

Introduction

Very recently, SiPMs were introduced in TD diffuse optics [START_REF] Dolgoshein | Status report on silicon photomultiplier development and its applications[END_REF]. These promising detectors are inexpensive and rugged, they do not require complex front end circuits, and they feature wide active areas of few square millimeters, thus permitting to maximize the photon collection. Essentially, they bring together the advantages of vacuum tubes and microelectronics detectors. After initial positive results in proof-of-concept studies [START_REF] Mora | Towards next-generation timedomain diffuse optics for extreme depth penetration and sensitivity[END_REF][START_REF] Mora | Fast silicon photomultiplier improves signal harvesting and reduces complexity in time-domain diffuse optics[END_REF], SiPM started to be integrated into compact and easy-to-operate detection modules [START_REF] Martinenghi | Time-resolved single-photon detection module based on silicon photomultiplier: A novel building block for time-correlated measurement systems[END_REF]. Thanks to the possibility to integrate into a single silicon chip both detectors and ancillary electronic circuits, SiPMs represent the most attractive technology for the next-generation of DOT systems [START_REF] Pifferi | New frontiers in time-domain diffuse optics, a review[END_REF]. Therefore, the demonstration of SiPM suitability also for TD DOT could open the way to novel immediate perspectives thanks to the possibility to conceive low-cost and reliable multichannel systems. Even if SiPMs have been already validated for diffuse optics in general, their suitability also for DOT is not straightforward. Indeed, as highlighted in [START_REF] Martinenghi | Spectrally Resolved Single-Photon Timing of Silicon Photomultipliers for Time-Domain Diffuse Spectroscopy[END_REF], the present generation of SiPMs employed as time-resolved single-photon detectors is affected by some limitations. The most critical one is the presence of a slow exponential decaying tail (known as "diffusion tail") detailed in section 3.1.1. Additionally, this long tail limits the dynamic range of the detector response function to about 2 orders of magnitude [START_REF] Martinenghi | Spectrally Resolved Single-Photon Timing of Silicon Photomultipliers for Time-Domain Diffuse Spectroscopy[END_REF], thus possibly affecting DOT reconstructions due to the reduced range suitable for data analysis [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF]. As reported in [START_REF] Martinenghi | Spectrally Resolved Single-Photon Timing of Silicon Photomultipliers for Time-Domain Diffuse Spectroscopy[END_REF], SiPMs are also characterized by high noise background (up to few hundreds counts per second), which limits the maximum number of photons per second that can be detected to avoid pile-up effects in TD acquisitions using time-correlated single photon counting systems [START_REF] Becker | The bh TCSPC Handbook[END_REF], and also the minimum photon bunch which can be extracted out of the noise.

In collaboration with the Politecnico di Milano, we explored the feasibility to perform timeresolved DOT with SiPMs and we demonstrate their performances to properly reconstruct, detect, localize and quantify the composition of scattering biological tissues in depth with a reflectance geometry.

For this study, we use the five totally absorbing objects defined in section 3.2.2.3 as absorption perturbations and we work with the same setup as the previous one with a source-detector distance fixed at 30 mm. First we remind and present the experimental setup and the signal analysis technique; then we analyze simulations and phantom measurements and we compare with results obtained with the well-known SPAD detectors. Finally we discuss the results and the perspectives.

Material and methods

Experimental setup

We used the same experimental setup as in section 3.2 except that we replace the detectors by SiPM detectors and we remove the timer in the schematized instrumental setup reported in Figure 3.7. So, this means that the light exiting from each collection fiber was then imaged directly onto a SiPM detector by a set of two lenses (achromatic doublets lens, focal length: 30 mm) providing a 1× magnification. Each module embeds a 1 mm 2 active area SiPM (C30742-11-050-T1, Excelitas Technologies, Canada) and all the required electronic circuits for sensing and amplifying the avalanche signal and yielding the pulse synchronous with the avalanche. For further details about the SiPM module see [START_REF] Martinenghi | Time-resolved single-photon detection module based on silicon photomultiplier: A novel building block for time-correlated measurement systems[END_REF]. Each SiPM module was then coupled to a TCSPC board (SPC-130, Becker & Hickl GmbH, Germany) to provide the start signal for timestamp generation. The stop pulse was given to both TCSPC boards by the "synchronization" output of the laser.

For the phantom, we used as previously the five totally absorbing inclusions of different volumes as perturbation objects. Their properties are reported in Table 3.2. The background medium was the same as previously with an absorption coefficient of 0.07 cm -1 and a reduced scattering one of 12 cm -1 at 820 nm. Scans of all the inclusions in the x-y plane were done at different depth (where "depth" refers to the distance from the surface of the phantom to the center of the inclusion) ranging from 10 mm up to 20 mm at step of 5 mm. Only for inclusion "size 3" also 25 and 30 mm depth scans were acquired in order to test the depth sensitivity of the setup. The depth was automatically changed using the motorized translation stage when a scan on the x-y plane was finished. The distance between injection and collection fibers was set to 15 mm and 30 mm depending on the experiment. In this section, we present only 30 mm source-detector distance. The results on the 15 mm interfiber distance probe are shown in Appendix C.

Preprocessing and analysis

The acquired or simulated measurements were preprocessed with the same approach as the SPAD study (paragraph 3.2.2.4) without the step concerning the gated mode. To simulate SiPM measurements, we took the same configuration as the experiments and we follow the simulation process defined in section 2.4.1 by time-convoluting computed Green's functions for each simulated phantom and the experimental IRF of the SiPM. For the reconstruction process, we apply the same method as the previous study; i.e. the standard code. However the space discretization was performed on a regular 80 × 65 × 45 mm 3 mesh grid of the medium because we probed object up to 30 mm depth.

Then the 3D reconstructions (coming both from measurements and simulations) was analyzed to extract spatial and absorption quantification information reported in Table 3.5. To get the values of spatial localization (x, y and depth), the centroids of each reconstructed inclusions were calculated on a delimited region defined by all values above a threshold of 50% of the recovered absorption variation ∆µ max a (where ∆µ max a is the variation between the reconstructed maximal absolute absorption value of the absorbing object and the one of the background). For quantitative analysis of these EBV inclusions, we use ∆µ vol a which is calculated by taking the average over a 1 cm 3 sphere at the expected depth of the inclusion. The expected values of the recovered x and y position are x = 0 and y = 0 for all depths. The quantification measures were given as relatives errors of ∆µ vol a in percentage in order to evaluate the accuracy of the retrieved values.

SiPM results

The results are presented into three parts: the simulation results which takes into account the experimental IRF of the SiPM, the phantom results and additional results on a comparison with the SPAD performances.

Simulation

Figure 3.15(a) shows the reconstructed maps (in the z-y and x-y planes, first and second row respectively) for the inclusion "size 3" at depths ranging from 10 mm up to 30 mm (columns).

The black circle reported in each map represents the volume of 1 cm 3 in which the absorption perturbation given by the inclusion is supposed to be. It is clearly noteworthy that in the z-y plane the position of the absorption perturbation is properly localized in depth till 25 mm, while for larger depths it is reconstructed in a position shallower than the real one. On the other hand, in the x-y plane there is a proper localization of the inclusion at all depths. For both z-y and Chapter 3. Detection and quantification of deep absorption perturbations using TR DOT 111

x-y slices at almost all depths, the reconstructed volume is larger than the expected 1 cm 3 and, more precisely, it progressively increases at higher depths. With the aim to study the sensitivity to the amount of absorption perturbation, Figure 3.16(a) reports the reconstructed maps (at a given depth: 15 mm) obtained using the whole set of totally absorbing objects. Regardless the volume, the position of the inclusion in both planes (z-y and x-y) is precisely reconstructed (for values see Table 3.5). In addition we clearly see that as the volume of the inclusion increases, the retrieved absorption perturbation ∆µ a increases while the reconstructed volume does not, in agreement with the model by Martelli et al [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 2: experimental implementation[END_REF][START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF]. However the retrieved value of the absorption perturbation is lower than predicted by Martelli et al. From simulations we can conclude that SiPMs are suitable for the tomographic than in simulations but this is most probably due to instabilities of the experimental system (e.g. laser instability, aging of the phantom) that were not taken into account in simulations. In any case, it is worth noting that the current setup is able to reconstruct an equivalent absorbing perturbation of 0.15 cm -1 (for 1 cm 3 volume) which is buried down to 25 mm in depth. Those values are typical for biomedical applications (e.g. optical mammography or brain imaging) thus making the use of SiPMs a good solution in medical imaging.

To ascertain the sensitivity on the amount of absorption perturbation, Figure 3.16(b) shows the z-y and x-y maps for the five inclusions at the same depth (15 mm). It is worth noting that all inclusions are properly reconstructed in term of depth and x-y position. Concerning the quantification of the absorption perturbation for experimental measurements, it is similar to that obtained from simulations for both the absolute value and for the trend showing a lower value of ∆µ a upon an increase of the depth.

As can be seen in Table 3.5, the depth localization is quite accurate for all depths with an error always smaller than 2.5 mm/4 mm (up to 20 mm/for 25 mm depth respectively) that corresponds to a maximum error of 17%. Most of time the depth at which the inclusion is reconstructed from both experiments and simulations are very close, although they may differ from the expected position but generally by less than 2 mm excepted at large depth (Figure 3.17). The retrieved coordinates of the x position of the inclusion recovered from the experiments seems to systematically differ by about 1 to 2 mm from the expected position. Maybe this is due to some systematic errors in the alignment of the inclusion in the experimental measurements.

On the other hand, the y localization presents a slighter error ( 1 mm) for both the simulation and the experimental data. To conclude, the error in the localization is small enough not to affect significantly the capability of the system to retrieve the position of an absorbing perturbation in a real application.

The inclusion "size 3" set at a depth of 30 mm produced a too low perturbation to be able to properly localize it not only in depth, but also in the x and y directions. For this reason, we can state that the maximum depth of the inclusion allowing proper tomographic reconstruction is 25 mm. Concerning the quantification capability, measurements and simulations are usually in good agreement (except for inclusion "size 1" and "size 2" at shallow depths). In any case, they are far from the expected value (reported in literature in [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 2: experimental implementation[END_REF][START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF]) and the error in quantification is higher upon an increased inclusion depth. Since this trend and the high discrepancy is present also in simulations, we suppose that it is not possible to ascribe it to an experimental problem but more to an issue connected with the model used for reconstruction or the IRF of the detector.

Comparison results with SPAD

To better evaluate the potential of the SiPM, we made a comparison with SPAD on the detection and quantification performances. Figure 3.17 presents quantified plots of the recovered depth of the small inclusion "size 3" in different configurations. By using SPAD in simulation, the gated mode enables to better detect at 30 mm depth with short source-detector distance. For 30 mm source-detector distance, both SPAD modes give comparable detection performances for small object at 10 to 30 mm depth (Figure 3.17(a)). Now, with the SiPM, Figure 3.17(b) shows some difficulties of detection starting from 20 mm depth for the experimental case at the short source-detector distance and starting from 25 mm depth for the three other cases.

If we compare both simulations using SPAD and SiPM at the large source-detector distance, we get the same behavior up to object at 20 mm depth then we see a faster decrease of the recovered depth for the SiPM.

Quantification

The same behavior for the quantification is observed for FG SPAD or SiPM results with ρ = 15 or 30 mm in simulation and in experiment (Figure 3.18). The same difficulty to quantify with the SiPM is observed when ρ = 15 mm (for the 3D reconstructions see Appendix C).

As already mentioned in the previous study with the SPAD, it is difficult to obtain accurate quantification for small, deeply embedded object of high ∆µ a . Both results using FG SPAD and

SiPM for inclusions at 10 mm depth show the limitation in the quantification ability although we can quantitatively characterize each inclusion: the plots behave similarly to the expected black curve and are linear up to ∆µ a = 15 cm -1 in all cases. The source-detector distance doesn't affect the quantification performances for object at 10 mm depth. These observations are confirmed by the simulations which give the same plots.

Discussion and conclusions

In this work silicon photomultipliers are proposed as potentially revolutionary detectors for time-domain diffuse optical tomography thanks to their unique features, bringing together the advantages of PMTs and SPADs. In particular, our aim was to validate their performance on heterogeneous phantoms to evaluate the impact of their inherent disadvantages: a slow tail and low dynamic range in the temporal response to single-photons, a high background noise and unknown stability for long measurement times in TCSPC applications. We performed a study based both on simulations (taking into account the true detector temporal response)

and experiments in order to disentangle non-idealities of the employed system from the actual detector behavior.

We demonstrated from both points of view that SiPMs clearly allow to detect perturbations of ∆µ a = 0.15 cm -1 in a volume of 1 cm 3 down to 25 mm depth into a scattering medium featuring µ a = 0.07 cm -1 and µ ′ s = 12 cm -1 using a 30 mm interfiber distance with good precision in lateral localization (better than 3 mm), while depth localization is excellent up to 20 mm (better than 2 mm). Additionally, it is worth pointing out that with a perturbation placed at a depth of 20 mm, a lateral localization better than 2.5 mm and a depth localization better than 2 mm can be achieved with a large range of perturbation values, ranging from ∆µ a = 0.056 cm -1 to ∆µ a = 0.94 cm -1 . By contrast, the accuracy in quantification of the absolute value of the perturbation is much worst, with a strong underestimation of the retrieved ∆µ a upon increasing depth. Still, this effect is not ascribed to the SiPM since it is replicated in simulations also for the FG SPAD.

The problem of the absolute quantification in diffuse optics measurements is well-known in the literature [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF] and a dedicated effort is needed to properly address this issue. Results with the SPAD (section 3.2.3) show that this issue is due to the physics of diffuse optics and the retrieval algorithm and not to the specific experimental implementation. By comparing with the FG SPAD performances, we notice similar detection performances up to 20 mm depth. Then the gated mode of the SPAD which increases the dynamic range permits to detect deeper inclusions with more accuracy as shown in Figure 3.17 and in [START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF]. Unfortunately, the present technology of SiPM does not allow the fast-gated operation. Concerning quantification performances at a given depth less than 20 mm depth, the SiPM and FG SPAD provide the same results.

The demonstration of SiPM suitability for TD DOT opens now the way to interesting perspectives. In the last chapter, we address the issue of quantification of multimaterial phantoms up to 20 mm depth with multispectral approach. So, we decided to use the SiPM at 30 mm source-detector distance because it has the same performances as the SPAD and it is potentially the future detector for TR DOT.

Chapter 4 MW TR DOT: material decomposition and concentration determination

For the final chapter, the whole multi-wavelengths TR reconstruction algorithm is applied on multimaterial phantom experiments at different depths. The purpose is to demonstrate the advantages of MW TR DOT to decompose and determine the chromophore concentration distribution inside biological tissues. Some current limits have been observed and discussed such as the degradation of the quantification with depth and the dependence with the initialization of optical properties in the reconstruction processing. The contribution of the spatially constrained method on the decomposition and the accuracy is confirmed.

In order to prepare the multispectral experiments, a multimaterial phantom study was carried out. Different absorbers and support phantoms (so called containers) were tested. Then a series of six gradual bimaterial inclusions were used as perturbed object and imaged at several wavelengths to retrieve the composition. We first paid attention to the absorption reconstruction at each wavelength and then to the material decomposition thanks to different combination of wavelengths.

Multimaterial phantom study

Optical phantoms are tissue-simulating objects used to mimic light propagation in biological tissues. They are mainly used for testing the performances of optical systems, calibration of systems, inter-lab comparisons and reference measurements [START_REF] Spinelli | Determination of reference values for optical properties of liquid phantoms based on Intralipid and India ink[END_REF]. In order to obtain accurate quantitative information, the manufacturing protocol has to be reliable, stable and well characterize with reference optical or non-optical methods.

In DOT studies, a variety of type of phantoms is used liquid, solid, gel, and mostly containing heterogeneities to investigate the system's ability to detect, delineate and quantify these inclusions.

Currently, the only standard phantoms are liquid phantoms based on Intralipid and inks. They are easy to use, adopted by most of the laboratories and a recent multi-center study demonstrates the high level of reliability to accurately characterize the optical properties of diffusive liquid phantoms at NIR wavelengths with an uncertainty of about 2% [START_REF] Spinelli | Determination of reference values for optical properties of liquid phantoms based on Intralipid and India ink[END_REF]. Nevertheless, optical phantoms containing heterogeneities are difficult to obtain and each lab has its own protocol [START_REF] Greening | Characterization of thin poly (dimethylsiloxane)-based tissue-simulating phantoms with tunable reduced scattering and absorption coefficients at visible and nearinfrared wavelengths[END_REF][START_REF] Pogue | Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry[END_REF]. One of the solutions to add inhomogeneities in liquid phantom is to use either a container filled with a liquid with chosen optical properties or a solid inclusion. But some complications can appear due to the effect of the container walls or the mismatch of the refraction index. These kinds of effects can be minimized by adapting the reference medium for example filling the inclusion with the same solution as the background medium.

In our framework of material decomposition, the heterogeneity must include several absorbers with proper and stable absorption spectra. These two challenges are studied below: first the existing absorbing substances, their optical spectra and the effect of the solvant on them, and then various possibilities to hold the inclusion.

Absorbers

Spectra of absorbers in water

The preparation of the multimaterial phantoms is critical and their optical properties have to be well-known a priori in order to compare with the reconstructed values. Beforehand, we have conducted a study to find an absorber with spectra similar to Hb and HbO 2 spectra and/or with well-distinguished spectra. All the spectra measurements were done with a spectrophotometer (Cary 300 UV-visible, Agilent). Figure 4.1 presents the spectral signatures of different colorants.

The blue, red, green and brown dyes are powders soluble in water from Gubra. The Luwsir dye is a NIR absorbing dye (dark green powder) from Luminochem. The cyan, black and magenta jet inks are basic inks for HP ink-jet printer. All the information about these additives are summarized in Table 4.1. Diluted in water, the spectra of the additives are stable in time except for the Luwsir dye which is not stable from one day to another; the absorption coefficient increases in the NIR range (light green curve Figure 4.1(b)). In this latter case, this would require to do control spectral measurements before and after the experiment. In Figure 4.1(a), the green dye has an interesting peak at 630 nm with a steep down slope and the blue dye has a specific "M" shape between 580 and 650 nm. The powder dyes can be mixed together and follow the law of mixtures. The jet inks have distinguished spectra but in the Vis range. The China black ink has a flat spectrum. We tested also food dyes sometimes used for manufacturing phantom [START_REF] Greening | Characterization of thin poly (dimethylsiloxane)-based tissue-simulating phantoms with tunable reduced scattering and absorption coefficients at visible and nearinfrared wavelengths[END_REF]. Their spectra are presented in Appendix D in presence of water or of ethanol. It is hard to find absorbers with an increasing slope after 700 nm. One good candidate would have been the Luwsir dye which follows very well the slope of the HbO 2 spectrum (680-750 nm) (Figure 4.1(b)) but it is stable only for a few hours. The most interesting optical contrasts are for lower wavelengths than the NIR range. So, it is probably more interesting to do acquisitions at wavelengths in the range of 550 nm to 700 nm by taking some cautions to remain in the validity of the diffusion theory and thus we should take an adequate value of the µ ′ s corresponding to these wavelengths. Figure 4.1(c) shows all the absorbers signatures with an isosbestic point at 600 nm. We can see on the left of this point that the steeper slope is from cyan jet ink and on the right it is the magenta jet ink but it is too vertical to be physically interesting; thus the other choice is the black jet ink.

Between 500 and 680 nm, the best choices to manufacture bi-absorber phantoms are cyan and black jet inks because their spectra are opposite variations with the wavelength and the fact that they are stable in solution for at least a week. Unfortunately, their spectra are far from those of the two hemoglobin forms (Figure 4.1(d)) but for our goal which is to validate the ability of TR DOT to decompose in chromophore it is not a constraint.

Effect of solvent on spectra

During the phantom study, we have tested inks diluted in ethanol. Actually dilution in ethanol is needed to mix it with polydimethylsiloxane (PDMS) used in the protocol of solid PDMS phantom described in the next section. By changing the solvent, significant modifications have been observed on spectra. In Figure 4.2(a), we diluted the jet ink in water and ethanol with the same proportions. Distortions and changes in absorption intensity are noticed on the three inks but especially for cyan jet ink. The significant effect of the solvent on spectra confirms the importance to characterize phantoms before and sometimes after the measurements. Moreover the PDMS itself strongly affects the absorption spectra of included dyes. To evaluate the effect of the PDMS stiff matrix, we made a mixture composed of transparent PDMS elastomer base with 1/10 of curing agent and some drops of a solution of jet ink and ethanol calculated to obtain µ a (600 nm) = 0.3 cm -1 . The mixture was poured into a disposable cuvette. As reference for the spectrophotometer, we took a cuvette filled only with the transparent mixture.

The spectroscopic measurements show that no similarity remains between ethanol-based solvent and PDMS-based solvent ( The conclusion of these observations is that every time a new solvent or material base to manufacture phantom is used, preliminary studies on the stability, reliably and optical characteristics must be done.

Phantoms for absorbent heterogeneities

After studying the absorber spectra, three different phantoms for inclusion are described below and tested in the multi-wavelength experiments in section 4.2. Figure 4.3 shows the pictures of the home-build support phantoms: resin container inclusion, PDMS inclusion and "balloon" inclusion.

A whole solid phantom (inclusion + background) seems to be preferable to avoid interface perturbations but it is not flexible, take time to built and it is difficult to characterize. In

Resin container inclusion

The first inclusion used was a resin container inclusion which previously existed in the lab. This inclusion was made of resin with added titanium dioxide (T iO 2 ) powder and black China ink to get µ ′ s = 12 cm -1 and µ a = 0.1 cm -1 at 750 nm. The resin object was machined to obtain a hollow resin 8 mm diameter and 8 mm length cylinder (inner dimensions) with thickness of 2 mm closed with two disks (waterproofed with Teflon) at the two ends. The thickness was large because of the fragility of resin. This kind of inclusion takes a long time to built and is difficult to characterize but it is stable in time and can be reused at will. It's also quite flexible about the content but there are some risks of bubble or leak which can soil the surrounding medium.

The thickness of the walls is an important problem because it can hide optical characteristics of the inner liquid phantom.

PDMS inclusion

Because a solid phantom is stable in time and easy to handle, we investigate the production of a series of PDMS-based phantoms. PDMS phantoms are permanent matrix compositions and are commonly used in diffuse spectroscopy as homogeneous medium that are suitable for routine calibration and for testing established systems [START_REF] Greening | Characterization of thin poly (dimethylsiloxane)-based tissue-simulating phantoms with tunable reduced scattering and absorption coefficients at visible and nearinfrared wavelengths[END_REF][START_REF] Roig | Caractérisation de tissus cutanés par spectroscopie bimodale : Réflectance Diffuse et Raman[END_REF].

To construct PDMS inclusions, we followed the mixture protocol developed in our laboratory by Dr. B. Roig [START_REF] Roig | Caractérisation de tissus cutanés par spectroscopie bimodale : Réflectance Diffuse et Raman[END_REF] and we used a plate of 48 holes (Falcon) (9.5 mm diameter and around 15 mm height each hole)(see Figure E.1(c) Appendix E) to obtain small 9.5×10 mm cylinder inclusions.

In parallel, the mixtures were also poured in plates of 6 holes (VWR) (35 mm diameter and 20 mm height each hole) in order to characterize them for more details Appendix E. Briefly to prepare these kind of inclusion, we first prepare the absorbing agent composed of ink and ethanol (strong concentration in order to avoid to put too much ethanol in the PDMS causing instabilities) and the scattering agent by mixing the curing agent and titanium dioxide (T iO 2 , 14027 Titanium(IV) oxide extra pure, Riedel-de-Haën). To determine the concentration of T iO 2 giving the wished diffusion coefficient, here µ ′ s = 13.5 cm -1 (600 nm), we use the theoretical formula of Dr. B. Montcel [START_REF] Montcel | Tomographie optique diffuse résolue en temps: Applications fonctionnelles en neurosciences[END_REF]. So, we added 520.48 mg of T iO 2 in 30 mL of curing agent and mixed with ultrasounds. Then we poured a volume of transparent PDMS elastomer base 1/10 of curing agent and we mixed by hand. The appropriate amount of absorbing solution should be added to the mixture. Before curing, the final mixture was mold in the plates with holes using a syringe and it must be degassed in vacuum for at least 90 min to remove air bubbles. To accelerate the curing, we placed the plate inside an oven at 70 ➦C during 2 h. The manufacture protocol is developed in details in the manuscript of Dr. B. Roig [START_REF] Roig | Caractérisation de tissus cutanés par spectroscopie bimodale : Réflectance Diffuse et Raman[END_REF].

The advantages of this base material are that scattering and absorption coefficients are tunable, the protocol is optimal to obtain the expected µ ′ s and the correct proportion of µ a but absorption accuracy is hard to obtain. It takes 2-3 days to manufacture a whole series of 6 inclusions but it is easy to embed them in the liquid phantom and they are reusable. Furthermore it was demonstrated above that the PDMS stiff matrix modifies the absorption spectra. The main manufacturing problems are the estimation of µ a and to have a homogeneous medium without T iO 2 sediments.

"balloon" inclusion

The third type of inclusion is called "balloon" inclusion. To manufacture "balloon" inclusions we cut the finger of a laboratory disposable glove (made of natural latex and powder free, Ansell Ltd, color clear) as support phantom. We attached the glove to a metal perforated plate and then the liquid phantom was blown inside though the hole. This inflating procedure allows to have a thicker support for phantom. We filled the support with 1 mL liquid phantom composed of distilled water, Intralipid and ink with a syringe of 10 mL. Then we made a knot with a white thread and fixed it to a white-painted wire in order to avoid any interference of the wire during the motion of the inclusion at different positions. The diameter of these inclusions was Correction for the effect of the support can be performed by filling the inclusion with the same solution as the background medium and using this as the "homogeneous" reference measurements [START_REF] Pogue | Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry[END_REF].

Multi-wavelengths experiments on a series of multimaterial phantoms

Thanks to these preliminary studies where we have chosen the two absorbers -cyan and black jet inks -and presented three different sort of inclusions, we can go to the next step the multiwavelength experimental TR-DOT measurements.

The experimental measurements were carried out at the Politecnico di Milano (LaserLab EU-ROPE funding) with a significant help of Dr L. Di Sieno, Dr A. Dalla Mora, Dr A. Farina and Prof. A. Pifferi. Different multimaterial phantoms with gradual percentage of cyan and black jet inks were tested with 3 different support inclusions based on resin holder, PDMS and "balloon". The aims of this campaign were to test new phantoms for multispectral measurements, evaluate our 3D material decomposition reconstruction in experimental conditions and assess the ability to quantify at different depths.

Instrumentation and phantoms

Instrumental setup

We used the same setup as in section 3.3 except that we change the laser source by a supercontinuum laser (SuperK EXTREME, NKT photonics, 10 ps pulse width) and add a fast wavelength selected system, an acousto-optical tunable filter (AOTF, NKT Photonics). The switching of the wavelength is fully automatized and controlled by a measurements software.

As in the previous setup, the backscattered light coming from the phantom was collected by two optical fibers at a distance of 30 mm from the excitation fiber (see Photo 4.4(a)). Each collection fibers were connected to two SiPMs coupled to a TCSPC system. Concerning the computation consideration, we used the Indirect method (See section 2.

Phantoms

The homogeneous liquid phantom was contained in a glass tank (a fish jar, volume 8.5 L) to avoid possible electrostatic problem which can happen with a plastic tank. It is composed of a water-based solution made of Intralipid and cyan jet ink (µ ′ s (600 nm) = 13.5 cm -1 and µ a (600 nm) = 0.095 cm -1 ) as exhibited in Picture 4.4(a). In order to assess the quantitation performances and to draw gradual quantified plots of absorbing heterogeneities, we realized five gradual bi-absorber inclusion compositions and one inclusion with the same solution as the background phantom. To have the optimal separation between the two colorants cyan and black jet inks, we decided to work in the optical window 550 to 660 nm and defined the proportions to get an isosbestic point of approximately 0.3 cm -1 at around 600 nm. Since we want to use the visible region, we used the "VIS channel" of the AOTF which covers a region from nearly 500 nm to 680 nm. The spectra of the gradual range of liquid phantom are represented in Figure 4.4(b). We started by 100% of cyan ink (µ a (600 nm) = 0.284 cm -1 ) up to 100% of black ink (µ a (600 nm) = 0.297 cm -1 ) with a decrease or increase step of 25% of each ink. The background phantom corresponds to around 33% of cyan ink and 0% of black ink (measures checked with the spectrophotometer).

The gradual absorbing heterogeneities were held up by the three different type of support inclusions:

1. The resin container inclusion of 8 mm diameter and 8 mm length was filled with the six bi-absorber liquid phantoms.

2. Six small cylinder gradual PDMS inclusions (9.5×10 mm) and six control PDMS phantoms were prepared following the protocol summarized in section 4.1.2. They consisted of T iO 2 and the same proportion of cyan and black jet inks as the liquid phantoms but we obtained flatter spectra as showed in We embedded within the homogeneous medium one after the other the 3 types of support inclusions with the 6 gradual compositions at 3 depths (10-15-20 mm). The depth is defined as the distance from the surface to the center of the inclusion.

Results

We present the results in three steps: first we compare the three support phantoms then we show and analyze the µ a reconstructions and finally the concentration reconstructions of the best type of inclusion.

Comparison of the phantoms

Beforehand, we compare the contrasts of the measurements of the three types of inclusion and evaluate if they are exploitable for the 3D DOT reconstructions. Figure 4.5 compares the contrasts of the acquisitions done above the position of the support inclusion at 10 mm depth containing the background medium (gray plots) and the 100% cyan liquid phantom (light blue plots) for the three sorts of inclusion and for each wavelength. The shortest wavelength 550 nm enables to probe deeper because the contrast can be visible up to around 4 ns before drowning in the noise in the case of the "balloon" inclusions. For the longest wavelength 660 nm, we can distinguish the contrast only up to no more than 2 ns.

Resin container inclusion

With the resin container phantom (left column), we observe a significant contrast of the container itself which decreases with the wavelength (gray plots), instead of a constant and very low contrast which would have been optimal. It seems to be caused by the thickness and the composition of the container (black China ink, slightly decreasing spectrum, see Figure 4.1(a)).

At 550 nm, we observe a big perturbation more than 60%. The effect of the container masks the spectral information of the liquid in the inclusion. The light blue contrasts (Figure 4.5) also decrease with the wavelength while we expected an increase of the contrast since cyan ink absorbs more at 600 nm than at 550 nm. We see a small separation between the two plots which would allow to detect the object in depth but too small to expect to quantify. Actually, we have attempted to reconstruct in 3D the distribution of the two absorbers in each inclusion (not shown). And as expected we can detect the inclusion but it is difficult to distinguish each inclusion and to quantify them. This resin container inclusion was not well adapted for this experiment (the chosen ink in the resin was not the same as in the background liquid phantom) and the thickness was too large but it is one of the limit of resin matter.

PDMS inclusion

Concerning the PDMS inclusion phantoms, we have shown above that the PDMS matrix affects the spectra of the included absorbers. So, the selected colorants were not adapted for our study because the spectra change in presence of PDMS, they are flatter and the difference between both are reduced. Consequently, it is difficult for the reconstruction to make a difference between the different compositions of the inclusions. By increasing the wavelength, we observe a decrease of the contrast of the PDMS inclusion supposing to have similar optical properties as the background medium. From 600 to 660 nm, the contrasts are close to 0 which means that the PDMS inclusion does not perturb measurements. Now if we look at the contrasts with the 100% cyan inclusion, less contrast is observed which reduces the chance to quantify something (Figure 4.5). We have tried to obtain the 3D reconstructions and we saw that we can detect at the 3 depths the inclusion but the quantification is unable to separate the inclusions (results not shown). With the chosen inks, the PDMS phantom gives hard exploitable contrasts and wrong results because of the non adapted spectra (too flat and low spectra) and the spectral difference with the surrounding medium. Maybe with other inks with more specific signatures in PDMS, these kind of phantoms can be used in DOT. The texture is interesting to mimic the skin and the PDMS is well used in spectroscopy. So, to built a whole phantom in PDMS can be interesting to have a probe coupling more realistic. The manufacturing protocol has still to be improved concerning the absorption calibration and compatible specific absorbing substances have to be found.

"balloon" inclusion

The contrast of the "balloon" inclusion itself (i.e. filled with background medium) is almost the same for all the wavelength (about 9%). So, the thickness of the glove seems to be wavelength independent and it will only add an offset in the measurements. The contrast of the 100% cyan inclusion increases from 550 to 600 nm and then it is stable ( Different contrasts have been observed for the other inclusions and the effect of depth shifts the beginning of the contrast at longer time as expected. In these experiments, the difference of contrast of the "balloon" inclusions is interesting and the "balloon" inclusion seems to be the best choice to accomplish quantitative analysis. This support phantom has less effect and is constant for each wavelength. For the rest of this chapter, we will only present results obtained with these inclusions where the measurements are the most exploitable and demonstrate the potential for nice quantitative reconstruction results.

µ a reconstructions and analysis

Thanks to the Indirect method for the reconstruction, we are able to detail the µ a quantification before the material decomposition. less absorbing inclusion. We observe the same kind of reconstruction artefacts around the red spots which forms a spreading reverse "L" probably due to the shape of the scan area.

In the z-y plane (Figure 4.7), all the inclusions show up in the reconstruction at the expected depth of 10 mm. A reconstruction artefact below the inclusion appears in every images with about the same shape. As for the x-y slices, the same gradation of the reconstructed composition is obtained. To test how the reconstruction depends on the depth of the inclusion, the z-y results of the "balloon" inclusion composed of 50% cyan and 50% black jet inks are drawn in Figure 4.8 for the 3 depths (10-15-20 mm). At 20 mm, the correct circular shape is lost and the red spot is more spreading due to the reduced SNR and the difficulty to predict the succession of scattering events to get information at deep layers. The reconstructed absolute absorption coefficient decreases by a factor that increases with the inclusion depth but it doesn't depend on the wavelength.

For the analysis of the accuracy of the quantification of the reconstruction, Figure 4.9 shows the maximal ∆µ a and the averaged ∆µ a over a 1 cm 3 spherical volume of each inclusion (markers) at each wavelength (colors: light blue, dark blue, purple, red and orange for 550, 580, 600, 630 and 550, 580 and 600 nm. For all the inclusions at 600 nm, we obtain values close to the expected all around ∆µ a = 0.2 cm -1 . At 630 and 660 nm, the results are overestimated in the range of -0.05 to 0.15 cm -1 (Figure 4.9(c)). The circle and triangle markers corresponding to the 0% cyan and 25% cyan ink inclusions respectively are the less accurate results. One of the reasons can be that their fraction of cyan ink is below that of the background phantom and thus below the initialization absorption coefficient µ

a . If we look at the oblique crosses corresponding to the inclusion filled with the background phantom, we reconstruct a constant ∆µ max a of around 0.04 cm -1 . Concerning the linearity, we get a coefficient of determination R 2 close to 1 for 550 nm and then we have also a correct linearity for 580 and 660 nm.

The values of the reconstructed ∆µ a averaged over a 1 cm 3 volume are smaller than the maximal ones but we keep the linearity for 550, 580 and 600 nm and the background inclusion has still constant value (Figure 4.9(b)).

With the spatially constrained method (i.e. forcing the reconstruction on the expected position of a sphere of 1 cm 3 ) (Figure 4.10), the underestimated values are raised and the overestimated values are decreased to come closer to the expected values (black line) . We gain on linearity and

Material decomposition

In the previous figures (rows 6 and 7 of Figures 4.6, 4.7 and 4.8), we have drawn the 3D concentration maps of the two absorbers in presence by taking into account all the wavelengths; i.e. the 5 wavelengths. As expected, the percentage of cyan ink increases gradually and the percent of black ink decreases gradually but we get an important underestimation with around 55% of black ink instead of 100%. For the background inclusion, we detect a small heterogeneity.

All the wavelength combinations have been tested and three of them are shown in Figure 4.11(a). We observe a good correlation with the cyan and black fractions with an increase of the percentage of cyan ink and a decrease of the percentage of black ink. For these wavelength combinations at 0% cyan ink, the quantity of cyan ink is really higher than expected which affects and thus reduces the measured fraction of black ink. Two explanations of this bias in the quantification are the initialization at 33% of cyan ink (background medium) and the offset given by the "balloon" container which both prevent to get reconstructed absorptions of the inclusion below the background optical properties. We see that for cyan fraction the closest results to the expected values are with the 550-580-600 nm combination (cyan continuous plot with the crosses). Concerning the black fraction, this latter combination (black continuous plot with the crosses) gives the best retrieved percentages for the three left points in Figure 4.11(a). Except for the 100% cyan inclusion where the best wavelength combination is with the 5 wavelengths, we can say in general that the best material decomposition is obtained with the wavelength combination 550-580-600 nm. Indeed these wavelengths yield us the best reconstructed ∆µ a (see Figure 4.9(a)). Using this optimal wavelength combination and the reconstruction with the spatial constraint (Figure 4.11(b)), we obtain very good correlation for both colorants at 10 mm depth (dotted line with points) with the expected values. For both colorants, they are close to the expected values for 50%, 75% and 100% cyan. At 0% cyan, we retrieve a value corresponding to the initialization (33 % cyan) and not below as we want. By increasing the depth, the linearity is kept but with a smaller slope coefficient and thus less accuracy. In comparison with the results obtained with the standard code with the inclusion at 10 mm depth, the a priori spatial information of the inclusion improves significantly the accuracy of the quantification. For the black fraction, we get accurate results with a maximal error of only 20% instead of 40% with the standard code when the object is at 10 mm depth. method using this optimum wavelengths. The gradation in the concentration of both absorbers is much more marked than by combining the 5 wavelengths. By forcing spatially the reconstruction (Figure 4.12(b)), the separation is evident and the quantification is encouraging. The problem of no negative reconstructed ∆µ a is still observed for the inclusions filled with 0% or 25% of cyan ink where we expect a dark blue (below 33%) spot instead of a green one.

Discussion and conclusion

Different multispectral phantoms have been tested and the development of "balloon" inclusion done in collaboration with the Politecnico di Milano is a promising, flexible and useful tool to evaluate the quantification capacity of a DOT system. Using this kind of inclusion filled with a gradual range of bi-absorber scattering solutions, we were able to reconstruct in 3D the composition of the probed phantom and separate distinctly the absorbers.

In laboratory experiment, the charaterization of phantom is primordial because we show that changing the solvant can modify at least the absorption properties and the thickness of a container can hidden the content information. Preliminar check are always necessary and enable to get rigorous protocol to evaluate the performances of the system.

The problem of the initialization is highlighted in this study. One interesting possibility to reduce this issue is to use the data taken at one of the wavelengths (typically the isobestic point of Hb/HbO 2 at 800 nm for biomedical application) as reference instead of an external reference.

This multispectral study demonstrates again the interest of the spatially constrained method ("a priori" code) to achieve the best of quantitative accuracy.

Conclusions and perspectives

Conclusions

The aim of this work was to evaluate in its globality the detection and quantification performances of multispectral time-resolved diffuse optical tomography in reflectance geometry in order to give quantitative medical images of the composition of biological tissues. We first worked on the reconstruction algorithm and on the analysis method. Then in simulation and experimental conditions, the quantification ability was investigated on single absorbing objects and further on multi-chromophore inclusions using a multispectral approach.

To achieve this goal, we tested in simulation two different approaches of processing the multiwavelength acquisitions in the reconstruction algorithm and we decided to use the "Indirect method" which is the most flexible. With this method, we obtain a 3D reconstruction map for each wavelength and thanks to the law of mixtures the concentrations can be retrieved. We also introduced a spatially constrained method which forces the geometry of the reconstruction by knowing a priori the size and position of the heterogeneities. The interest of this approach has been already demonstrated as an important improvement in the quantification in FD but only in the case of cylindrical geometry with a single inclusion at 10 mm depth in [START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF][START_REF] Srinivasan | Improved quantification of small objects in near-infrared diffuse optical tomography[END_REF]. In this PhD work, we extended the advantage of the constrained method in TD with a series of absorbing objects and we highlighted that at different depths up to at least 20 mm depth in a reflectance geometry this method improves the quantification.

Concerning the quantification performances, we report several conclusions. With numerical investigations, we have shown that our reconstruction algorithm is operational, robust and reliable to detect and reconstruct accurately the lateral (x and y directions) position of an inclusion. Additionally, we demonstrated that the quantification is not limited by the photonic noise but by the reconstruction provided that we apply a weighted least squares including the photonic noise to solve the inverse problem. The reliability of our simulation system has been confirmed from the agreement between simulations and experiments.

Concerning the experiment, we first evaluated the effect of the shape of the instrumental response function on the quantification capacity and we demonstrated that the IRF deteriorates significantly the contrast between the reference and the perturbed measurements when the exponential decay constant is approximately smaller than 10 ns Concerning the instrumentation, SiPM are promising detectors in TD systems [START_REF] Mora | Towards next-generation timedomain diffuse optics for extreme depth penetration and sensitivity[END_REF] thanks to their wide active areas permitting to maximize the photon collection and their low cost. Our feasibility study using these detectors in TR DOT has shown that the detection performances in the lateral localization are the same as the SPAD with less than 2 mm of error but the depth localization is limited to 20 mm depth because of the low dynamic range. For the quantification performances, we obtained the same results as FG SPAD with a correct quantification for the less absorbing objects at 10 mm depth although the absorption coefficient is underestimated at deeper depths and higher absorptions. This demonstration of SiPM suitability for TD DOT opens the way to interesting perspectives [START_REF] Pifferi | New frontiers in time-domain diffuse optics, a review[END_REF].

With these cheap detectors, we explored the chromophore decomposition thanks to multispectral measurements on a series of multi-chromophore inclusions. The behaviors of gradation of the two chromophores were in agreement with what was expected. For high concentrations, we still got an underestimation of the reconstructed values. In this study, the initialization with the known properties of the background medium affects significantly the quantification capacity.

The use of a spatially constrained method has again demonstrated its interest to achieve a more accurate quantification.

As a summary about the quantification performances, we demonstrated that the MW TR DOT has the potential to give quantitative medical images. Some understanding of the physics in diffuse optics has still to be improved and some enhancements in the reconstruction algorithm have to be done. Furthermore, the multimodal imaging can be a good opportunity to DOT for future clinical.

Perspectives

This PhD work on quantification has highlighted some important conclusions and clarified some difficulties by pointing out some tracks on their origins. We discuss in this section some proposals of solutions to overcome and improve the situation.

Multi-inclusions

To go to more complex media, an experimental campaign was carried out very recently on multi-inclusions phantoms. In collaboration with the Politecnico di Milano, we investigated bi-inclusions phantoms at a distance of 25 mm from their centroids (14 mm from their borders).

Preliminary results are displayed in Appendix F. These first results demonstrate the possibility to separate two objects at 25 mm distance with a 30 mm source-detector probe (and source position step of 7.5 mm), to view the percentage of the material composition, to identify and to localize them.

Algorithm improvements

In order to improve the algorithm, we suggest some possible improvements not yet tested. Like the three-steps reconstruction algorithm investigated in FD in [START_REF] Srinivasan | Improved quantification of small objects in near-infrared diffuse optical tomography[END_REF], we can proceed by steps.

One approach would be a two-steps reconstruction by first localizing the region of interest using the original algorithm and then perform a region-based reconstruction using the spatial prior information obtained from the first step. This second step would implement the Dictionary matrix D where the prior geometry information of step 1 would be injected in the reconstruction process. Another proposition, that keeps only one step, is to adapt the geometry of the reconstruction according to the clinical applications. For example, for flap surgery, we could modify the Dictionary matrix D to carry out a reconstruction by multi-layers instead of each voxels.

In the case of the small perturbations, an option is the Rytov approach presented in section 2.1.2.1. It could potentially be applied in one of the last stages of the reconstruction algorithm just before the solving of the inverse problem instead of replacing the Born approximation. It should be explored in more detail.

Otherwise, in the current algorithm there are some critical points which need to be improved like the stop condition in the iterations (currently fixed at 10 iterations). In order to have a more generic software, it should be better to find an objective criterion to stop the process when the convergence of the best estimation is reached. Another important point is the problem of the reference measurements and the initialization of optical parameters which is currently fixed from known external measurements. A possible solution should be to use a spectral reference (developed below).

Spectral reference

For in vivo and biomedical applications, one of the main problems of the perturbation approach in case of absolute measurement is the need of a reference medium. It is a critical problem.

Up to now, the reference measurements is carried out on a homogeneous phantom (with known optical properties) with a similar geometry as the perturbed medium. This method is restrictive.

With multispectral measurements, we are thinking about using a measurement acquired at one of the wavelengths as the reference measurement. In addition to being adapted to biomedical applications, another advantage of this method is that time is saved and the optical probe keeps the same position and coupling. The difficulties of this approach are that the IRFs must be the same for each wavelength. So, the time offsets (physical effect of the wavelength changing), the efficiency of the detector at the wavelength and the possible differences in shape should be taken into account by acquiring beforehand the IRFs. With our "cross iterative Born approximation", we need to know the initial optical properties of the reference medium and we consider the reference medium as homogeneous. In the spectral reference approach, the reference medium does not fulfill these criteria but the possible solution is to make a deconvolution to retrieve the

A.1 Results

As in section 3.

B.3 Conclusion

The performances of non-gated and fast-gated SPAD was addressed in this Appendix. With the FG SPAD, we get more focused, less scrambled 3D reconstructions and it is possible to detect small object of 5 mm diameter at least up to 30 mm depth. Till 20 mm depth, the non-gated and fast-gated SPAD 3D reconstructions are the same with both source-detector distances. For short source-detector distance, the gated mode brings improvements in the detection of deep objects which is interesting and important to design small held-hand probe.

For quantitative analysis of these EBV inclusions, we analyze the relative errors of ∆µ vol a in percentage in order to evaluate the accuracy of the retrieved values (Table C.1). ∆µ vol a is calculated by taking the average over a 1 cm 3 sphere at the expected depth of the inclusion.

Table C.1 are reported the quantification errors in simulation and in experiment. At 10 mm depth, the quantification is the best with around 40% of relatives errors in average by taking all the inclusions but with only the three first inclusions we get a smaller errors of 25% in average.

For deeper positions, the quantification is worse with errors going up to 99%. spectra. They absorb in the visible range to around 690 nm. From 550 to 700 nm, the red food dye spectra decreases with a steep slope. The blue food dye has a peak at 630 nm in presence of water or of ethanol. For these kind of dye, no shape difference is observed between the water or ethanol solvent but they are not absorbing in the NIR range. The spectra are normalized at 600 nm to obtain 0.3 cm -1 .

D.2 Spectra of Lunir dye

The Lunir dye (reference: LUNIR5 from Luminochem) is comparable to Luwsir dye because they are NIR absorbing dye and are both a dark green powder but Lunir dye is soluble in ethanol contrary to Luwsir dye which is only soluble in water. In order to obtain a homogeneous mixture, we put the beaker in a water chamber and mix a few minutes with ultrasounds. This dye gives an interesting spectrum high absorbing and increasing in the NIR range. So, the Lunir dye can be an attractive absorber to built PDMS phantoms for DOT but a stability study in ethanol and PDMS base is first needed.

µ a (600 nm) = 0.3 cm -1 and ends with an inclusion containing 100% cyan ink with the same µ a (600 nm) = 0.3 cm -1 . The gradation for the intermediates is by steps of 25% of each ink. To construct PDMS inclusions, we followed the mixture protocol developed in our laboratory by Dr. B. Roig [START_REF] Roig | Caractérisation de tissus cutanés par spectroscopie bimodale : Réflectance Diffuse et Raman[END_REF]. To produce the series of inclusion, we first prepared the absorbing agent composed of ink and ethanol (strong concentration in order to avoid to put too much ethanol in the PDMS causing instabilities) and the scattering agent by mixing the curing agent and titanium oxide (T iO 2 , 14027 Titanium(IV) oxide extra pure, Riedel-de-Haën). To determine the concentration of T iO 2 giving the desired reduced scattering coefficient, here µ ′ s = 13.5 cm -1 (600 nm), we used the theoretical formula of Dr. B. Montcel [START_REF] Montcel | Tomographie optique diffuse résolue en temps: Applications fonctionnelles en neurosciences[END_REF]. So, we added 520.48 mg of T iO 2 in 30 mL of curing agent and mixed with ultrasounds. Then we poured a volume of transparent PDMS elastomer base (Sylgard 184 (Dow Corning); refraction index = 1.41 at 700 nm), 1/10 of curing agent and we mixed by hand. The appropriate amount of absorbing solution should be added to the mixture. Before curing, the final mixture was mold in the plates with holes using a syringe and it must be degassed in vacuum for at least 90 min to remove air bubbles. To accelerate the curing, we placed the plate inside an oven at 70 ➦C during 2 h. We used a plate of 48 holes (Falcon) (9. The inclusion 2 is more cyan and the inclusion 3 is more gray which is in agreement with the color in water.

E.2 Characterization

To check the optical properties, we tested two methods of characterization.

Spectrophotometer

It was demonstrated by Dr. B Roig [START_REF] Roig | Caractérisation de tissus cutanés par spectroscopie bimodale : Réflectance Diffuse et Raman[END_REF] that her protocol is reliable to obtain the expected µ ′ s . For the µ a , we verify the calculated proportion of ink with the transparent PDMS elastomer base before adding the scattering agent by measuring it with a spectrophotometer. To do that, we make a mixture composed of the transparent PDMS elastomer base with 1/10 of curing agent and some drops of a solution of each jet ink and ethanol calculated to obtain µ a (600 nm) = 0.3 cm -1 . The mixture is poured into a disposable cuvette ( signatures of black and cyan jet inks and thus leads to more difficulties to separate the two absorbers. The PDMS solid matrix seems to change the chemical structure affecting the absorption properties.

Time-resolved spectroscopy

With a time-resolved spectroscopy approach, we are potentially able to evaluate the reduced scattering and absorption coefficients. These characterization measurements were carried out by Nous avons également mis en place une méthode de contrainte spatiale en adaptant le code en incorporant a priori la position et la taille de la zone d'intérêt, dans notre cas l'inclusion. Nous explorons cette approche avec de l'information a priori afin d'améliorer la quantification.

Dans la quatrième section du chapitre, nous avons évalué la robustesse et la fiabilité de l'algorithme pour la reconstruction, la détection et la quantification des milieux diffusants. La première étude de simulation a permis de tester différentes adaptations du programme. Elle a démontré que le bruit de photons n'affecte pas les résultats de reconstruction grâce à la bonne prise en compte de celui-ci dans l'algorithme. Une détection précise au niveau des directions latérales est obtenue et la détection en profondeur est également bonne grâce à la pondération des voxels par le carré de la distance entre le voxel et la position source-détecteur introduite dans l'algorithme.

En simulation, la méthode de contrainte spatiale permet de retrouver la cartes de µ a initiales. La seconde étude de simulation compare les deux méthodes de l'approche multispectrale et c'est la méthode Indirecte qui est retenue en raison d'une légère meilleure quantification et de sa flexibilité pour combiner différemment les longueurs d'onde.

Conclusions et perspectives

Le but de ce travail de recherche était d'évaluer dans sa globalité les performances de détection et de quantification de la tomographie optique diffuse multispectrale et résolue en temps en géométrie de réflectance dans le but de donner des images médicales quantitatives de la composition des tissus biologiques. Nous avons d'abord travaillé sur l'algorithme de reconstruction et sur les méthodes d'analyse. Ensuite, en simulation et dans les conditions expérimentales, la capacité de quantification a été étudiée sur de simples objets absorbants et sur des inclusions multi-chromophores avec une approche multispectrale.

Les conclusions de ce travail sont que la DOT multispectrale résolue en temps a le potentiel pour produire des images médicales quantitatives. Nous sommes capables de quantifier de façon linéaire des faibles absorptions à 10 mm de profondeur et la décomposition en matériaux donne des résultats suivant une tendance en accord avec les proportions des chromophores. Les points faibles qui ont été identifiés sont d'une part la sous-estimation des valeurs quantifiées augmente avec la profondeur et d'autre part la forte influence de l'initialisation de l'algorithme avec les caractéristiques du milieu de référence qui peut biaiser la quantification.

En ce qui concerne l'apport de l'a priori dans l'algorithme de reconstruction, c'est un véritable succès pour l'amélioration de la quantification. L'imagerie multimodale peut être une bonne opportunité pour que la DOT puisse s'insérer dans le milieu clinique.

Plusieurs perspectives sont envisageables pour développer la MW TR DOT. La référence spectrale permettrait de s'affranchir d'une référence extérieure. Des améliorations du code de reconstruction sont à essayer telles qu'une approche plus systématique pour arrêter le reconstruction avec un critère de convergence et une approche de reconstruction en deux étapes avec d'abord la localisation puis la quantification.

Abstract

In medical imaging, the ability to accurately retrieve and quantify the composition of turbid media is challenging and would enable to diagnose some diseases or to better study physiological processes. Diffuse optical tomography (DOT) is an attractive medical imaging technique which permits to probe in depth using near-infrared light and to reconstruct in three dimensions the composition of biological tissues non-invasively, non-ionizing and with potentially high specificity. To obtain endogenous chromophore (oxy-and desoxy-hemoglobin) features in the depth of a highly scattering medium, a multiwavelength time domain optical setup combined to a three-dimensional reconstruction algorithm was developed. Experimental measurements were conducted in reflectance geometry by illuminating a perturbed medium (with a heterogeneity) with a picosecond laser and by collecting, for several wavelengths and multipositions, the backscattered light via two fibers connected to two dedicated detectors and coupled to a time-correlated single photon counting system. The data processing of these time-resolved measurements and those of a known reference medium was performed by supposing that the propagation of light is governed by the diffusion approximation and using a method based on Mellin-Laplace transform. Numerical and phantom experiments on series of objects similar to biological media demonstrate that this technique has the potential to give quantitative medical images. We have highlighted a correct quantification for the less absorbing objects at 10 mm depth while underestimation results at deeper depths and higher absorptions. Furthermore, the multimodal imaging brings improvements in quantification in depth and thus it can be a good opportunity to DOT for its future clinical applications.

Keywords: Diffuse optics, signal processing, reconstruction, time-resolved detection, multispectral, quantification Résumé La possibilité de déterminer exactement et de quantifier la composition des milieux biologiques est un défi pour l'imagerie médicale qui permettrait de diagnostiquer certaines maladies ou de mieux étudier les processus physiologiques. La tomographie optique diffuse (DOT) est une technique d'imagerie attrayante qui permet de façon non invasive, non-ionisante et potentiellement avec une grande spécificité de sonder les milieux en profondeur en utilisant la lumière dans le proche infrarouge et de reconstruire en trois dimensions leur composition. Pour obtenir les caractéristiques des chromophores endogènes (oxy-et désoxy-hémoglobine) enfouis dans un milieu très diffusant, une instrumentation optique dans le domaine temporel et multi-longueurs d'onde combinée à un algorithme de reconstruction en trois dimensions a été développée. Des mesures expérimentales ont été menées en géométrie de réflectance en éclairant avec un laser picoseconde un milieu perturbé (contenant une hétérogénéité) et en récupérant, pour plusieurs longueurs d'onde et multi-positions, la lumière rétrodiffusée via deux fibres optiques connectées à deux détecteurs dédiés et couplés à un système de comptage de photon unique. Le traitement des données de ces mesures résolues en temps a été réalisé en supposant que la propagation de la lumière est gouvernée par l'approximation de la diffusion et en utilisant une méthode basée sur la transformée de Mellin-Laplace. Des simulations et des expériences sur une gamme de milieux imitant les milieux biologiques ont démontré que cette technique médicale a le potentiel pour donner des images médicales quantitatives. Nous avons montré que l'on peut déterminer correctement la composition d'objets à 10 mm de profondeur absorbant faiblement. Pour de plus fortes profondeurs et des absorptions plus élevées, la valeur du coefficient d'absorption ou de la concentration est sous-estimée. En outre, l'imagerie multimodale apporte des améliorations dans l'exactitude de la quantification en profondeur et donc peut être une bonne opportunité pour les futures applications cliniques de la DOT.

Mots clés : Optique diffuse, traitement du signal, reconstruction d'images, détection résolue en temps, multispectrale, quantification

ε

  i (λ)C i (1.3)

Figure 1 . 2 :

 12 Figure 1.2: The three regimes for DOT: continuous-wave (CW), frequency-domain (FD) and time-domain (TD) where detected light intensity over time is illustrated on (1), (2) and (3) respectively. On the left, schematic distribution of photons in a sample in the reflection and transmission geometries. As a rule of thumb, the mean light penetration depth in the reflection geometry is of order ρ/2 (Extracted from[START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF]).

Figure 1 . 3 :

 13 Figure 1.3: TCSPC measurement principle [10]. The original waveform is reconstructed by detecting single photons, which are randomly distributed during the laser pulse time duration, and by recording photon counts by means of an histogram graph.

  labs on DOT reconstructions. Other research groups have developed their own homemade software reconstruction. A non-exhaustive list of the active groups on DOT reconstructions are in the USA the team of Dr. D. Boas at the Massachusetts General Hospital and Harvard Medical School, the team of D. A. H. Hielsher of Columbia University with their VITO system, in Europe the group of Dr. H. Wabnitz at Physikalisch-Technische Bundesanstalt (Germany), the group of Dr. H. Dehghani at University of Birmingham, the group of Dr. J. Hebden at UCL, the group of Dr. A. Liebert at Institute of Biocybernetics and Biomedical Engineering (Poland) and Dr. L. Hervé at CEA-LETI (France) and in Asia Dr. Feng Gao at Tianjin University (China), Dr. Koichi Shimizu at Hokkaido University (Japan).

Figure 1 . 4 :

 14 Figure 1.4: (a) FD DOT reconstructed images coming from[START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF] of an absorbing cylindrical object (8 mm diameter) at z = 0 mm (x = -22 mm and y = 0 mm). The images represent slices from z = -50 to 50 mm with a step of 10 mm. The x and y axis correspond to the horizontal and vertical directions respectively. The top row doesn't include a priori information on the position and shape of the object contrary to the bottom row (figures from[START_REF] Dehghani | Threedimensional optical tomography: resolution in small-object imaging[END_REF]). (b) and (c) FD DOT quantitative results of the average reconstructed absorption coefficient plotted as a function of the times factor to the background of perturbed objects of 10 mm and 20 mm diameter (b) after using the 3 steps reconstruction algorithm and (c) after using region-based reconstruction with a priori information on the position and size of the object (Data extracted from[START_REF] Srinivasan | Improved quantification of small objects in near-infrared diffuse optical tomography[END_REF]).

( a )Figure 1 . 5 :

 a15 Figure 1.5: Reconstructed images from clinical applications. (a) NIR mammography images (Figure from [97]). (b) neonatal brain hemorrhage images (Figure extracted from [49]). (c) vascularization foot images and vessel localizations (Figure from [60]).

Figure 1 . 6 :

 16 Figure 1.6: Reconstructed images from brain activity on different language areas using highdensity DOT and fMRI. The last row shows the overlap between these two techniques (Figure extracted from [38]).

  diagnosis. Multimodal imaging has been first investigated by the combination of positron emission tomography (PET) and X-ray computed tomography (CT) or magnetic resonance imaging (MRI) and ultrasound imaging (US) ... Recently, diffuse optical imaging systems were implanted by research groups in multimodal imaging platforms such as MRI/DOT, CT/DOT, US/DOT, photoacoustic tomography (PAT)/DOT and other possible combinations provided complementary information (Chapter 10 of [65]).

Figure 1 . 7 :

 17 Figure 1.7: (a) Hybrid combined US/NIR system and a handheld probe with NIR optical fibers at the periphery of the US array developed in the University of Connecticut. (b) a. US image of indeterminate breast lesion in 37-year-old woman. b. Total hemoglobin (HbT) concentration maps represent into slices from 0.3 cm underneath the skin surface with 0.5 cm spacing (in absolute). The maximal HbT is 35.3 µmol/L. So it is a benign lesion. (c) a. US image of indeterminate right breast lesion in 65-year-old woman. b. HbT concentration maps represent into slices from 0.5 cm underneath the skin surface with 0.5 cm spacing. The maximal HbT is 112.3 µmol/L. So it is a malignant lesion. Biopsy revealed an invasive carcinoma. (Figure extracted from [133]).

Chapter 2 Development

 2 of a multi-wavelengths (MW) time-resolved (TR) DOT three-dimensional reconstruction algorithm In this chapter, we explain the multi-wavelengths time-resolved DOT (MW TR DOT) reconstruction algorithm developed at our laboratory. In the first part, forward model and inverse

Figure 2 .

 2 Figure 2.1 shows the diagram of the reconstruction process in case of monochromatic acquisitions. The initial optical properties of a reference medium, the spectra of the expected material in the medium, the geometry of the source-detector positions and the geometry of the mesh are injected in the direct model to generate the computed measurements. Then these model data Y model are compared with the acquired data Y exp . If the difference between Y model and Y exp converges, or if a predefined number of iterations is reached, the reconstruction process is stopped. Otherwise, the differences between measurements (Y exp ) and predictions (Y model ) are used to update the optical properties through a step called "inverse problem". At the end of the reconstruction process, we obtain optical properties maps which are compatible with the measurements.

Figure 2 . 1 :

 21 Figure 2.1: Reconstruction algorithm. Y model and Y exp correspond to the model and acquired measurements respectively. µ (k) a is the absorption coefficient distribution at iteration k. N loop is the stop condition. X is the solution obtained from the inverse problem.

Figure 2 . 2 :

 22 Figure 2.2: (a) and (b) Plots of the temporal Green's functions for a source-detector couple (ρ = 15 mm) generated though a heterogeneous medium composed of various absorbing objects (1 cm 3 sphere) at 10 and 20 mm depth respectively (log scale). (c) and (d) The corresponding contrasts between the reference Green's functions (black curve) (homogeneous medium) and the perturbed Green's functions with inclusion at 10 and 20 mm depth respectively.

Figures 2. 2

 2 Figures 2.2(a) and (b) show the generated TR Green's functions of several heterogeneous media

  Figures 2.2(a) and (b) show the generated TR Green's functions of several heterogeneous media

  wavelengths and at several source/detector positions. Time-resolved signal is acquired from the backscattered light in a reflectance geometry. These data are used in the reconstruction algorithm (Figure 2.1).

Figure 2 . 3 :

 23 Figure 2.3: Illustration of the time convolution of the time-resolved measurements. Here the IRF is the IRF of fast-gated single photon avalanche diode (FG SPAD).

Figure 2 . 4 :

 24 Figure 2.4: Normalized instrumental response function (IRF) in log scale of a fast-gated single photon avalanche diode (FG SPAD) with the representation of the three main characteristics of an IRF: decay (slope in log-scale), dynamic range and full-width at half maximum (FWHM). This IRF has a decay equal to exp(-9.4 × t), a dynamic range of 6 decades and a FWHM of 258 ps.

  Figure 2.5(a) illustrates the relation between the arrival times of detected photons and their travel area in the medium.

Figure 2 . 5 :

 25 Figure 2.5: (a) Schema of the time-resolved perturbation approach. (b), (c), (d) and (e) timeresolved signals of a fast-gated single photon avalanche diode (FG SPAD) with a 15 mm sourcedetector distance probe on the surface above the position of a totally black inclusion (5 mm diameter equivalent to ∆µ a = 0.15 cm -1 ) at 5, 10, 15 and 20 mm depth. The red curves correspond to the reference TR measurements through the reference medium A, M Aλ sd (t). The blue curves correspond to the perturbed TR measurements through the perturbed medium B, M Bλ sd (t) (log scale). The time step is 4.1 ps/channel.

1 .

 1 As already mentioned, a study has been done by Dr. A. Puszka (section 2.1.3 and 2.1.4 of her PhD manuscript[START_REF] Puszka | Diffuse optical tomography: a time-resolved approach for reflectance measurements at short source-detector separation[END_REF]) on the comparison of different datatypes (Moments and Mellin-Laplace transforms). The conclusion is that the MLT offers a versatile formalism which enables to easily switch between different datatypes. The MLT is a time-filter which enables to time-window the signal reducing the time dimension and increasing the sensitivity for late time where photons collected have traveled longer distances. The MLT has two parameters which give the shape of the windows: n corresponding to the index of the number of orders N (n : 0 → N ) and p corresponding to the number of extracted transforms by nanosecond. As in Dr. A. Puszka's PhD manuscript we fix p = 3 ns -1 and we will go to N = 20 which is reasonable for extracting information at late time (i.e. in depth) and reasonable computation times.

Figure 2 . 6 .

 26 photons and, therefore, to improve the quality of reconstruction of deep layers of diffusive media. So, it is important to use high orders in order to probe in depth. By looking at the contrast plots in Figure2.6, the contrast is important also at late time (Figure2.6(f)) and thus high orders and a large enough p (p = 3 ns -1 ) enable a comparable contrast with the temporal contrast (Figure2.6(h)). So the MLT has kept the main information and permits to work on

Figure 2 . 6 :

 26 Figure 2.6: (a) and (b) MLT windows with p = 1 ns -1 and p = 3 ns -1 respectively (N = 20).(c) Reference and perturbed TR measurements using FG SPAD with inclusion (1 cm 3 sphere of ∆µ a = 0.15 cm -1 ) at 10 mm depth and (f) its contrast. (d) and (e) Transformation of the TR measurements through the MLT and (g) and (h) their corresponding contrasts for p = 1 ns -1and p = 3 ns -1 respectively.

Figure 2 . 7 :

 27 Figure 2.7: Simulation process.

Figure 2 .

 2 Figure 2.8: (a) Geometry of the simulated experiment and phantom. (b) and (c) time-resolved simulated measurements of the reference and the perturbed medium above the inclusion with and without Poisson noise at ρ = 15 mm and 30 mm respectively. (d) and (e) contrasts at the 3 depths (10-15-20 mm) with and without noise for ρ = 15 mm and 30 mm respectively.

  Figure 2.8(b) and (c) show the reference and perturbed TR simulated measurements without and with Poisson noise for 15 mm and 30 mm source-detector distancerespectively. We can observe that the FWHM of the signal is larger for ρ = 30 mm which is explained by the fact that the photons have to travel a longer distance and thus the time-offlight is longer and the contrast appears later. The presence of the inclusion at 10 mm depth modifies the perturbed signal which is smaller than the reference measurements because photons have been absorbed by the absorbing object.

  photons arrive sooner and the probability of detection is higher.

Figure 2 .

 2 Figure 2.10: 3D reconstructions of simulation study 1 using standard 2 code. On the top, reconstruction with inclusion ∆µ a = 0.15 cm -1 at 15 mm depth and on the bottom at 20 mm depth.

Figure 2 . 11 :

 211 Figure 2.11: Overlaid profiles of the 3D reconstructions for the 5 codes in all the configurations (ρ = 15 or 30 mm, without or with inclusion, at 10-15-20 mm depth, without or with noises).

Figure 2 .

 2 12(b) shows the reference and perturbed measurements above the position of the inclusion 11 with λ = 750 nm. In this section, we tested the two multispectral methods: Method 1 is the Indirect algorithm and Method 2 is the Direct algorithm. And we applied the standard code with the MLT parameters: p = 3 ns -1 and n = 20. The visualization of the 3D reconstruction images (Figure 2.13) displays the Hb and HbO 2 concentration decomposition. We also draw the HbT concentration by summing Hb and HbO 2 and the percent of oxygen saturation ratio (SO 2 ) obtained by calculation with Eq. 1.8. Quantification analysis. To estimate the performances of quantification of the different multispectral methods, we analyze the reconstructed data with 4 different quantification techniques (Figure 2.14): in concentration we extract the maximal reconstructed concentrations of Hb and

Figure 2 .

 2 Figure 2.12: (a) Schematic of the simulated phantom and scan configuration with ρ = 15 mm. (c) Absorption spectra of the 11 simulated inclusions.(b) time-resolved measurements. (d) contrasts at each wavelength above the inclusion 11 at 10 mm depth.

Figure 2 .

 2 Figure 2.12(d) presents the contrast of inclusion 11 at different wavelengths. In agreement with the spectrum of inclusion 11 (see cyan curve 2.12(c)), the contrast decreases from 750 nm to

Figure 2 . 13 :

 213 Figure 2.13: 3D reconstructions of Hb, HbO 2 , HbT and SO 2 for inclusion 6 (50%) and inclusion 11 (100%) with Method 1: Indirect method (λ by λ). Left column slices z-y and right column slices x-y. Black circles and rectangles are the true inclusion positions.

Figure 2 .

 2 Figure 2.13 shows the 3D reconstructions of Hb and HbO 2 and also HbT and SO 2 for inclusion 6 and inclusion 11 with the same colorbar. The reconstruction method used for these reconstructions is the Indirect method (Method 1). The positions of the inclusion are correct and has the expected shape except for HbO 2 images where the center is located a little too deep in z. Small spreading problems in depth (z axis) can be observed. The retrieved concentrations are increasing from inclusion 6 to 11 and the background values are close to the expected values (C Hb = 25 µM and C HbO 2 = 32 µM ). As expected, the percent of SO 2 is smaller than in the surrounding medium. Concerning the accuracy, the HbT concentration is underestimated and

Figure 2 . 14 :

 214 Figure 2.14: Reconstructed quantification plots with different quantification techniques: (a) maximal reconstructed concentrations Hb and HbO 2 , molar quantities (b) by integrating over the expected volume of the inclusion (a priori ), (c) by taking the integral over an area defined with a threshold of 30% of δC max and (d) by integrating over the expected volume of the inclusion + extension of 3 mm. Method 1: Indirect (λ by λ) and Method 2: Direct algorithms.

Figure 2 .

 2 Figure 2.14(a), we see that with Method 1 the maximal reconstructed Hb concentration slightly increases and is accurate up to inclusion 5. With the technique using the integral over the a priori volume of the inclusion, the best method is Method 1 where the oxyhemoglobin quantities are correct. The other techniques (integral over a support volume defined as 30% of ∆C max and

  photon can be safely detected. The active area is few tens of micrometers in diameter. One real advantage of SPADs compared to other devices is the capability to switch between ON and OFF states very quickly (in less than 1 ns). As we shall see, this can be used to increase the dynamic range of the TCSPC signals (Figure3.1).Actually, the Politecnico di Milano conceived their own time gating mode although such performance is not available in commercial devices. These detectors are called fast-gated (FG) SPAD[START_REF] Boso | Fast-gating of single-photon avalanche diodes with 200 ps transitions and 30 ps timing jitter[END_REF][START_REF] Tosi | Fast-gated single-photon counting technique widens dynamic range and speeds up acquisition time in time-resolved measurements[END_REF]. A proof-of-principle done by A. Puszka et al[START_REF] Puszka | Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes[END_REF] confirms the potential of FG SPADs to image deep absorbing inclusions in TR DOT. They offer the shortest rise time, a good sensitivity in the NIR range, a very thin IRF and a reasonably flat gate and increase imaged depth by using fast-gating. Planar silicon SPADs offer high quantum efficiency in the visible spectrum and in the first part of the NIR range (about 60% at 550 nm, 20% at 800 nm and still 4% at 950 nm).

Figure 3 . 1 :

 31 Figure 3.1: Preprocessing on SPAD measurements. (a) The 5 recorded measurements at 5 consecutive delays. (b) Then each acquisition is multiplied by the corresponding attenuation according to the delay. (c) The full TPSF reconstructed for gated (black curve) and non-gated (dashed blue curve) SPAD (cut curve for the non-gated curve because of negative values). (d)The signal to noise ratio (SNR) of gated and non gated SPAD (board 2)(log scale).

Figure 3 . 2 : 8 *

 328 Figure 3.2: Theoretical and experimental IRFs in (a) log scale and in (b) linear scale. α is the exponential decay constant of theoretical IRFs.

Figure 3 . 3 :

 33 Figure 3.3: Effect of IRF down slope on the contrast obtained by different 1 cm 3 spherical absorbing objects (∆µ a = 0.056, 0.087, 0.15, 0.37, 0.94 cm -1 ) at different depths (10-15-20-25-30 mm) with fixed background medium at µ a = 0.07 cm -1 . The source-detector distance is 15 mm.

3

 3

  and 3.4, in the ideal case with an "ideal" IRF (black curves), the deeper is the object, the later the contrast starts and the slower is the increase of the contrast. It is mainly due to the travel of the photons which is longer to meet a deep inclusion and to come back which is caused by the scattering phenomenon. The contrast increases logarithmically as a function of the ∆µ a (exact equation: Contrast% = 92.2 + 46.16 × log 10 (∆µ a )). For a same absorption perturbation, the change of the background absorption has no influence on the contrast (Figure3.4). Compared to the other contrast plots, we observe a very small offset in time which is due to the small FWHM of 4.1 ps allocated on the theoretical IRFs.

Figure 3 . 4 :

 34 Figure 3.4: Effect of the IRF decay on the contrast by changing the background medium µ a = 0-0.07-0.14-0.2 cm -1 . The source-detector distance is 15 mm.

Figure 3 . 5 :

 35 Figure 3.5: Effect of the 5 experimental IRFs after background subtraction on the contrast obtained by different 1 cm 3 spherical absorbing objects (∆µ a = 0.056, 0.087, 0.15, 0.37, 0.94 cm -1 ) at different depths (10-15-20-25-30 mm) with fixed background medium at µ a = 0.07 cm -1 .

3 ,

 3 we present two experimental campaigns done at the Politecnico di Milano with the funding of a LASERLAB-EUROPE grant (grant agreement no. 654148, European Union's Horizon 2020 research and innovation programme) and with the important investment of Dr L. Di Sieno, Dr A. Dalla Mora and Dr. A. Farina. Both experimental investigations used the same experimental protocol except that the first one was done with non-gated and fast-gated SPAD detectors and the second one with SiPM detectors.

I

  wish to thank Dr. L. Di Sieno, Dr. A. Dalla Mora, Dr. A. Farina and Prof. A. Pifferi of the Dipartimento di Fisica of the Politecnico di Milano for theirs investments in the experimental part conducted in Milan with their detectors dedicated for time-resolved acquisition and their instrumentation setup.

Figure 3 .

 3 6 (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 = -40 mm). Photo 3.6 (b) and (c) show the two experimental geometric configurations.

Figure 3 . 6 :

 36 Figure 3.6: (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 mm. y = 0 mm). The crosses are the detection positions of the couple of detector optical fibers at 15 mm (blue) and 30 mm (green) distance from the source. Each black dot is separated by 7.5 mm in both directions. (b) and (c) Photos of the two probe configurations. (b) Photo with the source fiber (yellow arrow) and two detectors (blue arrows) at 15 mm interfiber distance fixed in the black PVC plank and the tank is filled with the liquid phantom composed of Intralipid and black ink. Here the phantom is ready for acquisition. (c) Photo with interfiber distance 30 mm between the source position (yellow arrow) and the two detectors (green arrows).The phantom is in preparation because the tank is half-empty to change the inclusion.

  7) 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 through an optical fiber (core diameter: 200 µm; numerical aperture (NA): 0.22; 2.45 m long; Lightech srl, Italy), as reported in Figure 3.7. Photons reemitted from the sample were collected in two different positions by means of two fibers (core diameter: 1 mm; NA: 0.37; 2.45

Figure 3 . 7 :

 37 Figure 3.7: Instrumental setup with the phantom containing an absorbing object (detailed in section 3.2.2.3). Schematic experimental setup of the two studies. A light pulse generated with a tens of picosecond pulsed laser at a defined wavelength goes through a variable optical attenuator (VOA) before illuminating a phantom with an absorbing inclusion at a given depth. Then the backscattered light is collected with two detectors coupled to two TCSPC boards connected to the computer. An electronic synchronization signal at the output of the laser is sent to the TCSPC boards. The green spots correspond to the scan positions of the probe (one source and two detectors).

Figure 3 . 8 :

 38 Figure 3.8: Photo of the 5 totally black cylinders. From left to right, inclusions with 3.2, 4, 5, 6.8 and 8.6 mm diameter. The equivalent ∆µ expected a

2 . 3

 23 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 ∆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 ∆µ a . For high values of N L, the slope must be evaluated for the effective x value. A lower value of b -accompanied by a low non-linearity N L -indicates still a linear trend but with a reduced slope on the retrieved ∆µ a . Obviously, for b → 0 and N L = 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 , thus for a null absorption perturbation. The closeness of b to 100%, accompanied by a low and coefficients can be used as another measure of accuracy.

  and 3.10 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.

Figure 3 . 9 :

 39 Figure 3.9: Numerical and experimental 3D reconstruction absorption maps with non-gated SPADs at 30 mm 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 = 10 mm). From top to bottom, the 4 inclusion sizes are shown (see Table3.2 for equivalent ∆µ a ). The black circle corresponds to a 1 cm 3 sphere and is centered on the expected position of the inclusion.

Figure 3 .

 3 10 shows the effect of the depth on quantification. The reconstructed absorption of a given inclusion ("size 2" is shown here) decreases with the depth (10-15-20 mm). The experimental and simulation 3D maps are comparable and follow the same trend. In both Figures (Figures 3.9 and 3.10), 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 (see Annexe B).

Figure 3 . 10 :

 310 Figure 3.10: Numerical and experimental 3D reconstruction absorption maps with non-gated SPADs at 30 mm 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).

Figure 3 . 11 :

 311 Figure 3.11: Overview of quantified plots of ∆µ vol a , 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 10, 15 or 20 mm 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.

Figure 3 . 12 :

 312 Figure 3.12: Relative errors in percent (%) of ∆µ vol a between (a) simulation and reality (Sim/Th) and (b) experiment and reality (Exp/Th) calculated with Eq. 2.43 for each reconstruction in each configuration (ρ = 15 or 30 mm and Non Gated (NG) or Gated (G) SPAD). Colors encode depth while the filling gated or not and the color intensity the ρ.

Figure 3 . 13 :

 313 Figure 3.13: 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 10, 15 or 20 mm 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.

Figure 3 .

 3 Figure 3.14: [Constrained method] Relative errors in percent (%) of ∆µ a between (a) simulation and reality (Sim/Th) and (b) experiment and reality (Exp/Th) calculated with Eq. 2.43 for each reconstruction in each configuration (ρ = 15 or 30 mm and Non Gated (NG) or Gated (G) SPAD). Colors encode depth while the filling Gated or Not Gated and the color intensity the ρ.

Figures 3 .

 3 Figures 3.15 and 3.16 show the z-y and x-y slices of the 3D reconstructions in simulation and experiment using SiPMs. It permits to compare side by side the similarities and discrepancies.

  (a) Simulation: slices of 3D reconstructions. (b) Experiment: slices of 3D reconstructions.

Figure 3 .

 3 Figure 3.15: z-y and x-y slices of 3D reconstructions of a perturbation equivalent to inclusion "size 3" at different depths with 30 mm interfiber distance (a) in simulation and (b) in experiment .

  (a) Simulation: sections of 3D reconstructions. (b) Experiment: sections of 3D reconstructions.

Figure 3 .

 3 Figure 3.16: z-y and x-y slices of 3D reconstructions of the 5 perturbation objects at 15 mm depth with 30 mm interfiber distance (a) in simulation and (b) in experiment.

Detection

Figure 3 . 17 :

 317 Figure 3.17: Quantified plots on the detection performances (recovered depth in function of the true depth) with an inclusion of 5 mm diameter (inclusion "size 3") embed at 10 to 30 mm depth. (a) in simulation using non-gated (NG) and fast-gated (FG) SPAD with ρ = 15 or 30 mm. (b) using SiPM in simulation (simu) or in experiment (exp) with ρ = 15 or 30 mm (source-detector distance). The black curves correspond to the expected curves.

Figure 3 . 18 :

 318 Figure 3.18: Quantified plots on the quantification performances (recovered ∆µ vol a in function of the expected ∆µ expected a ) with the 4 first inclusions embed at 10 mm depth. Comparison between results obtained from FG SPAD or SiPM with ρ = 15 or 30 mm (source-detector distance) (a) in simulation or (b) in experiment.

Figure 4 . 1 :

 41 Figure 4.1: Spectra of different absorbers in presence of water (a) with a common point at 800 nm and (c) with a common point at 600 nm. (b) Spectra of different dyes compare with the ones of Hb and HbO 2 with isosbestic point at 800 nm. (d) Spectra of black and cyan jet inks compare with the ones of Hb and HbO 2 with isosbestic point switched at 600 nm. (µ a = 0.3 cm -1 is expected at 600 nm).

  Figure 4.2(b)). The PDMS matrix flattens the optical signatures of black and cyan jet inks and thus leads to more difficulties to separate the two absorbers. The PDMS solid matrix seems to change the chemical structure affecting the absorption properties as demonstrated with the spectra in the Vis-NIR range in Figure 4.2(b).

Figure 4 . 2 :

 42 Figure 4.2: Spectra of jet inks for printer (a) in presence of water and ethanol (b) and also PDMS where µ a = 0.3 cm -1 is expected at 600 nm.

  other imaging techniques (magnetic resonance imaging and ultrasound imaging), agar-based phantoms have been used for decades. However to have more flexibility and fast manufacturing series of phantoms, we decided to work with surrounding liquid phantom composed of water, Intralipid and ink and solid or container phantoms filled with liquid solution as inclusions.

Figure 4 . 3 :

 43 Figure 4.3: Pictures of phantom supports for inclusion: (a) resin container inclusion filled with liquid phantom, (b) solid inclusion built with PDMS and (c) "balloon" inclusion filled with liquid phantom.

11. 5

 5 mm on average. Because of the thin thickness of the inflating glove, we can perceive the color of the phantom by transparency (Figure F.1(a) of Appendix F).

Figure 4 . 4 :

 44 Figure 4.4: (a) Picture of a part of the setup with the probe fixed at a PVC plank (ρ = 30 mm) and the phantom. (b) Liquid phantom spectra for each inclusion measured via the spectrometer before adding scattering agent and the five selected wavelengths (550-580-600-630-660 nm) for the scanning (red lines).

  2.1.1) for practical reasons. It is more flexible to use this technique because we can test different combination of wavelengths. Moreover we did preliminary tests in Chapter 2 with the two methods and we have concluded that to fully analyze the quantification potential the Indirect method is the optimal choice. In fact, we have found weighting problems with the Direct methods. About the experimental time and the computation consuming time, a whole scan with 35 source positions and two detectors, 5 wavelengths, 3 depths and 5 repetitions gives us 120 min of acquisition time then for the preprocessing 1.3 s for one depth/one wavelength and 5 min for the reconstruction of one scan at one wavelength and one depth.

  Figure 4.2(b) and explained in section 4.1.1.2. 3. Six "balloon" inclusions were manufactured by taking care to carry out the same protocol for each inclusion, injecting the same quantity of substance (1 mL) and checking their height and diameter; we measured 11.5 mm diameter on average with a standard deviation of 0.6 mm. Both types have been filled with the six gradual liquid phantoms (Figure 4.4(b)).

Figure 4 . 5 :

 45 Figure 4.5: Contrasts above each inclusion's type at each wavelength at 10 mm depth containing the same solution as the background medium ("bkg" gray line) and containing the solution 100% of cyan ink ("100% cyan" light blue line). At long time, the gray curves are hidden by the light blue curves.

Figure 4 . 5 )

 45 in agreement with the spectrum measured with the spectrophotometer before adding the Intralipid (Figure 4.4(b)).

Figures 4 .

 4 Figures 4.6and 4.7 display the 3D reconstructed maps of the "balloon" inclusions at 10 mm depth on two slices y-x and y-z respectively (z being the coordinate along the depth with the origin at the surface of the scattering medium). The five first rows correspond to the absorption maps for each wavelength and the two last rows which will be analyzed in the next section are the % of both inks maps deduced from the mixture law (Eq. 1.3) taking into account their respective absorption spectra.

Figure 4 . 6 :

 46 Figure 4.6: Slices y-x of 3D reconstruction maps of each "balloon" inclusion at each wavelength at 10 mm depth and the material decomposition in percent of cyan and black inks by taking into account the 5 wavelengths (rows 6-7). Same colorbar for each line.

Figure 4 . 7 :

 47 Figure 4.7: Slices z-y of 3D reconstruction maps of each "balloon" inclusion for each wavelength at 10 mm depth and material decomposition in percent of cyan and black inks by taking into account the 5 wavelengths. Same colorbar for each line.

Figure 4 . 8 :

 48 Figure 4.8: Slices z-y of 3D reconstruction maps of 50%-50% "balloon" inclusion for each wavelength at 3 different depths (10-15-20 mm) and material decomposition in percent of cyan and black inks by taking into account the 5 wavelengths. In this figure, the colorbars are defined independently for each image.

Figure 4 . 9 :

 49 Figure 4.9: Comparison of the reconstructed and expected values of ∆µ a for all the "balloon" inclusion (markers) and and various wavelengths (colors) at 10 mm depth. The linear fits of each wavelength are plotted (dotted lines and equation of the fit). (a) ∆µ a maximal. (b) ∆µ a averaged over a 1 cm 3 spherical volume at the expected depth; 10 mm. The shape of the symbols corresponds to the composition of the inclusions, their color to the wavelengths. (c) Relative errors on ∆µ max a . (d) Relative errors on ∆µ vol a .

Figure 4 .

 4 Figure 4.10: [constrained method] Comparison of the reconstructed and expected values of ∆µ max a for all the "balloon" inclusions (markers) and different wavelengths (colors) for (a) 10 mm (same as Fig. 4.9(a)) and (b) 15 mm inclusion depth when imposing in the reconstruction a spatial constraint for the position and size of the inclusion. The shape of the symbols corresponds to the composition of the inclusions, their color to the wavelengths. The linear fits of each wavelength are plotted (dotted lines). Relative errors on ∆µ max a at (c) 10 and (d) 15 mm depth.

Figure 4 . 11 :

 411 Figure 4.11: Comparison of the reconstructed and expected maximal cyan and black fractions (%) according to the "balloon" inclusion at (a) 10 mm depth with three different combination of wavelengths and (b) at all the depths with the constrained method with the optimal combination of wavelengths (550-580-600 nm).

Figure 4 .

 4 Figure 4.12: (a) Slices z-y of 3D material decomposition in percent of cyan and black inks of each "balloon" inclusion at 10 mm depth by taking into account wavelengths 550, 580 and 600 nm and (b) with constrained method.

Figure 4 .

 4 Figure 4.12 show the results of the 3D material reconstruction without and with the constrained

- 1 .

 1 So, to avoid important effects of the IRF, the decay constant has to be considered in the specifications of the choice of the detectors in the TR DOT chain. The fast-gated SPAD developed by the Politecnico di Milano are the best current detectors thanks to their steep decay and also to their high dynamic range with enables to detect deep objects. The IRF of the SiPM has a similar decay but it has a slow tail at long time and a low dynamic range. From experimental and simulation results, we got accurate and linear quantification of small inclusions embedded at 10 mm depth from 0 to 0.25 cm -1 in absolute absorption. For more absorbing objects, a deviation from linearity and a decrease of the accuracy is observed. The depth also affected the accuracy of the retrieved absorption coefficients or concentrations by underestimating the expected values (by a factor 2, 3 and 6 for z = 10, 15 and 20 mm respectively). The spatially constrained method completely solves the problem in simulation and gives reliable accurate results in experiment. For future biomedical applications, this approach is not unrealistic, since it can be implemented in coregistration with a different imaging modality yielding the size and location of brain activations or suspect lesions.

Figure A. 1 :

 1 Figure A.1: Theoretical IRFs in (a) with FWHM = 260 ps and in (b) with FWHM = 781 ps (log scale). α is the exponential decay constant of the theoretical IRFs.

1 . 2 . 1 ,

 121 we convolute the simulated Green's function of various phantoms composed of various 1 cm 3 absorption perturbations (∆µ a = 0.056, 0.087, 0.15, 0.37, 0.94 cm -1 ) at five different depths(10-15-20-25-30 mm) with the theoretical IRFs. The background medium was kept at 0.07 cm -1 in absorption and 12 cm -1 in scattering. The obtained contrasts between the reference and the perturbed TPSFs using the seven defined IRFs is represented in Figure A.2 for the FWHM of 260 ps and in Figure A.3 for the FWHM of 781 ps. In Figure A.2, we observe almost the same results as in Figure 3.3 with the same contrast values at long time. Only for the 10 mm deep inclusions we can see a short transient where we still see a increase of the contrast according to the increase of absorption perturbation and the effect of the depth which delays as expected the appearance of the contrast. So, between 4.1 ps and 260 ps, the FWHM doesn't modify the contrast that depends mostly on the slope of the IRF.The FWHM of 781 ps simulates the impact of a bad FWHM. Figure A.3 presents the results on the contrasts. Now, the short transient mentioned above where all the curves were superimposed extends to a longer duration such that it can be visible even for the 20 mm inclusion depth.It results from the large FWHM but no alteration at long time is obtained. The quantitative information in deep layer is conserved and it can be recovered from the reconstruction.

Figure A. 2 :

 2 Figure A.2: [FWHM = 260 ps] Effect of IRF down slope on the contrasts obtained by different 1 cm 3 spherical absorbing objects (∆µ a = 0.056, 0.087, 0.15, 0.37, 0.94 cm -1 ) at different depths (10-15-20-25-30 mm) with fixed background medium at µ a = 0.07 cm -1 .

Figure A. 3 :

 3 Figure A.3: [FWHM = 781 ps] Effect of IRF down slope on the contrasts obtained by different 1 cm 3 spherical absorbing objects (∆µ a = 0.056, 0.087, 0.15, 0.37, 0.94 cm -1 ) at different depths (10-15-20-25-30 mm) with fixed background medium at µ a = 0.07 cm -1 .

Figure B. 1 :

 1 Figure B.1: z-y and x-y slices of 3D reconstructions of a perturbation equivalent to inclusion "size 3" at different depths with 15 mm interfiber distance (a) using non-gated SPADs and (b) using fast-gated SPADs both in simulation. The black circle corresponds to a 1 cm 3 sphere and is centered on the expected position of the inclusion.

Figure B. 4

 4 Figure B.4: z-y and x-y slices of 3D reconstructions of the 5 perturbation objects at 10 mm depth with 30 mm interfiber distance (a) in simulation and (b) in experiment using the FG SPAD. The colorbar is common for each column.

  (a) Simulation: sections of 3D reconstructions. (b) Experiment: sections of 3D reconstructions.

Figure C. 2

 2 Figure C.2: z-y and x-y slices of 3D reconstructions of the 5 perturbation objects at 10 mm depth with 15 mm interfiber distance (a) in simulation and (b) in experiment using SiPM.

Figure D. 2 shows

 2 the spectrum of lunir dye diluted in ethanol which was done on the spectrophotometer.

Figure D. 2 :

 2 Figure D.2: Absorption spectrum of Lunir dye in ethanol.

Figure E. 1 :

 1 Figure E.1: Photos of the protocol to built PDMS phantoms. (a) The twelve mould PDMS phantoms for the characterization in two plates of 6 holes (VWR) (35 mm diameter and 20 mm height each hole). The first row is similar to the second row. (b) The six final PDMS inclusions used in the experimental campaign of Chapter 4. (c) Plate of 48 holes (Falcon) (9.5 mm diameter and around 15 mm height each hole) used to mould the inclusions. (d) Three cuvettes filled with transparent PDMS base, ink and ethanol for the absorption measurements with a spectrophotometer.

  5 mm diameter and around 15 mm height each hole)(Figure E.1(c)) to obtain small 9.5×10 mm cylinder inclusions. In parallel, samples for the optical characterization of the material were prepared by pouring the same mixture in plates of 6 holes (VWR) (35 mm diameter and 20 mm height each hole)(Figure E.1(a)).It takes 2-3 days to manufacture a whole series of six inclusions.

  Figure E.1(b) shows the six final PDMS inclusions used in the experimental campaign. We see a color difference by eyes.

  Figure E.1(d)). As reference for the spectrophotometer, we took a cuvette filled only with the transparent mixture (note "Ref" in Figure E.1(d)).

Figure E. 2

 2 Figure E.2(a) displays the spectroscopic measurements of the 100% cyan and 100% black inks represented by the dotted lines. We observe that these spectra are different compared to the spectra of the two absorbers in water (dashed lines). The PDMS matrix flattens the optical

Figure E. 2 :

 2 Figure E.2: Spectra of PDMS phantoms composed of cyan and black inks. (a) Absorption spectra of the expected values in water (dashed lines), the measured values in transparent PDMS base with a spectrophotometer (dotted lines) and the measured values in scattering PDMS with time-resolved spectroscopy (continuous lines). (b) Scattering spectra of the normalized theoretical values based on Prahl's works [76] (dotted line), the theoretical values based on Montcel's formula [77] (dashed line), and the measured values in scattering PDMS with timeresolved spectroscopy (continuous lines).

  Dr A. Farina from the group of the Politecnico di Milano in transmittance geometry with their instrumentation on the larger PDMS phantoms mould at the same time and thus with the same mixtures as the inclusions. They were made especially for the characterization (Figure E.1(a)). The method of characterization for retrieving absorption and scattering is based on analytical model under the diffusion approximation. They fit the convolution between the model and the IRF to the data using a Levenberg-Marquardt non linear fitting algorithm. They minimize the reduced χ 2 . Usually the fitting range is 80% of the peak on the rising edge and 1% of the peak on the falling edge. The results are shown in Figure E.2 by the continuous lines. For the scattering properties (Figure E.2(b)), the theoretical values of the plot of Montcel and of Prahl are similar from 630 to at least 830 nm. The formula of Montcel was used in his works only from 690 to 870 nm

Figure F. 3 :Chapitre 1 :

 31 Figure F.3: Profiles obtained from the black lines on the slices of Figure F.2 by taking into account the 3 wavelengths (550-580-600 nm). The left inclusion is localized at x = -20 mm and y = 0 mm (continuous plots). The right inclusion is localized at x = 5 mm and y = 0 mm (dashed plots). The blue and black plots correspond to the profiles of the percentage of cyan and black inks respectively. The red lines correspond to the depth position of each inclusion.

Chapitre 3 :

 3 Détection et quantification de perturbation d'absorption en profondeur en utilisant la DOT résolue en temps Avant de faire des campagnes de mesures expérimentales, nous avons évalué l'impact de l'instrumentation sur la quantification. Plusieurs réponses instrumentales (IRF) simulées et mesurées ont été testées en les convoluant à des fonctions de Green générées (solution de l'équation de diffusion lorsque la source est un Dirac c'est à dire la réponse idéale du milieu) de plusieurs milieux perturbés par une inclusion. Il a été conclu que la pente descendante définie par la constante de décroissance exponentielle de l'IRF est un facteur important pour réaliser de corrects mesures quantitatives et lorsque cette constante est inférieure 10 elle commence à avoir un effet destructeur sur le contraste. Le meilleur détecteur actuellement sur le marché est le SPAD déclenché qui est développé par le Politecnico di Milano et qui a une décroissance raide et aussi une plage dynamique élevée permettant de détecter des objets profonds. Ensuite, deux campagnes d'expériences ont été menées sur la détection et la quantification d'objets absorbants en profondeur en utilisant la DOT résolue en temps avec deux types de détecteurs: les SPADs et les SiPMs. Pour ces deux études, nous avons reproduit en simulation les expériences et des résultats similaires ont été obtenus. D'après les résultats expérimentaux et de simulations avec les SPADs, nous avons obtenu une quantification précise et linéaire pour les petites inclusions placées à 10 mm de profondeur de 0 à 0.25 cm -1 en absorption absolue. Pour des objets absorbants plus fortement, une déviation de la linéarité et une diminution de la précision sont observées. La profondeur affecte également la précision des coefficients d'absorption récupérés en sous-estimant les valeurs attendues (par un facteur 2, 3 et 6 pour z = 10, 15 et 20 mm respectivement). La méthode de contrainte spatiale résout complètement le problème pour la simulation et donne des résultats fiables et précis. Pour les applications biomédicales futures, cette approche n'est pas irréaliste, car elle peut être mise en oeuvre en la combinant avec une autre modalité d'imagerie donnant la taille et l'emplacement de l'hétérogénéité.En ce qui concerne la deuxième campagne, les SiPMs sont des détecteurs prometteurs dans les systèmes optiques du TD grâce à leurs larges zones actives permettant de maximiser la collecte de photons et à leur faible coût. Notre étude de faisabilité de l'utilisation de ces détecteurs pour la DOT résolue en temps a montré que les performances de détection et de la localisation latérale sont les mêmes que avec le SPAD avec moins de 2 mm d'erreur, mais la localisation en profondeur est limitée à 20 mm de profondeur en raison de la faible plage dynamique. Pour les performances de quantification, nous avons obtenu les mêmes résultats qu'avec les FG SPADs avec une bonne quantification pour les objets faiblement absorbants à 10 mm de profondeur, puis une sous-estimation à des profondeurs plus profondes et des absorptions plus élevées. Cette démonstration de l'aptitude des SiPMs à être utilisée pour la TD DOT ouvre la voie à des perspectives intéressantes.

Chapitre 4 :

 4 Tomographie optique diffuse multispectrale résolue en temps : décomposition sur une base de matériaux et détermination des concentrations Avec les SiPMs, nous avons exploré les méthodes de détermination de la composition d'inclusions multichromophores grâce à des mesures multi-longueurs d'onde sur une série de fantômes liquides. Tout d'abord, nous avons fait une étude sur des fantômes multi-matériaux et testé plusieurs techniques de réalisation d'inclusion. Des inclusions "balloons" faites avec le doigt d'un gant de laboratoire gonflé avec un fantôme liquide (eau, intralipide et encre) ont été définies comme les meilleurs modèles pour obtenir des inclusions multi-matériaux car l'épaisseur du gant est faible et conduit à des perturbations petites et constantes. L'importance de la caractérisation des fantômes a été mis en évidence en observant des variations de spectres d'absorbants en fonction du solvant utilisé comme l'eau et l'éthanol. Pour la campagne d'expériences, nous avons fabriqué une gamme de six inclusions "balloons" avec une fraction graduelle d'encre noire et d'encre cyan (encres pour imprimante HP). En décomposant sur une base de matériaux, les reconstructions en fonction des variations de proprotions d'encres cyan et noire montrent une tendance conforme à ce qui était attendu. Pour les fortes concentrations, nous retrouvons une sous-estimation des valeurs reconstruites comme dans les expériences du chapitre précédent. Nous avons également observé que l'initialisation par les propriétés optiques du fond empêche la reconstruction de valeurs en dessous de ce fond et donc entraine des biais dans la quantification. En revanche, la méthode de contrainte spatiale améliore fortement la quantification mais elle est aussi affectée par le problème de l'initialisation.

  is a software package which simulates light propagation in biological tissue based on the finite-element method (FEM) for FD and CW available online http://www.dartmouth.edu/ ~nir/nirfast/.

	This work was developed at Dartmouth College for FD NIR imaging. NIRFASTSlicer, a
	recent interative gui driven toolbox, enables analysis, visualization, segmentation of med-
	ical images, modeling and image reconstruction in optical imaging. Since a few years,
	the group of Dr. H. Dehghani organizes a workshop on NIRFAST jointly with the OSA
	Biomedical Optics Congress.

❼ TOAST++ (temporal optical absorption and scattering tomography)

[START_REF] Schweiger | The Toast++ software suite for forward and inverse modeling in optical tomography[END_REF]

is the opensource image reconstruction package software for time-resolved DOT created at the University College of London available online http://web4.cs.ucl.ac.uk/research/vis/ toast/demos_matlab.html. TOAST++ is the last update version of TOAST. This soft-

  It uses a CW approach and gives optical images of blood flow without compression of the breast. ART Advanced Research Technologies Inc. (Canada) had started in 2008 to sell SoftScan an optical breast imaging system TR spectroscopy but the company failed in 2013. Its non-invasive and painless approach has used TD optical imaging technology. It might have allowed clinicians to better locate and characterize breast tumors as benign or malignant and could have provided faster assessment of therapeutic effectiveness.

	Florida (USA) founded in 1993 but is not yet approved by USA Food and Drug Administration
	(FDA).
	CTLM (Computed Tomog-
	raphy Laser Mammography) is commercialized by Imaging Diagnostic System Inc. located in

Table 2 .

 2 2: Summary of x, y and z in mm extracted from the 3D reconstructions (a) without noise and (b) with noise added in the measurements and theirs standard deviations on the 10 noise draws in parentheses.

				x (stdx)	y (stdy)	z (stdz)
	noise depth	code	sd15	sd30	sd15	sd30	sd15	sd30
			stand 2	0.12	0.17	-0.16	-0.42	12.33	11.58
			stand 1	0.02	0.13	-0.05	-0.18	10.28	9.64
		10	R1 2	0.04	0.14	-0.24	-0.51	12.47	11.55
			R1 1	0.05	0.15	-0.09	-0.19	10.86	9.80
			a priori	0	0	0	0	10	10
			stand 2	0.07	0.39	-0.04	-0.19	17.29	17.22
			stand 1	0.02	0.15	0.08	0.06	14.52	15.03
	no	15	R1 2	0.07	0.34	-0.06	-0.24	18.02	17.61
			R1 1	0.04	0.16	0.06	0.03	15.22	15.52
			a priori	0	0	0	0	15	15
			stand 2	0.02	0.14	0.07	-0.03	21.56	21.10
			stand 1	0.01	-0.02	0.16	0.21	17.98	17.73
		20	R1 2	0.02	0.12	0.02	-0.07	22.56	21.85
			R1 1	0.04	-0.01	0.16	0.18	18.64	18.04
			a priori	0	0	0	0	20	20
					(a) without noise
			stand 2 0.10(0.00) 0.20(0.00) -0.16(0.00) -0.43(0.00) 12.25(0.01) 11.44(0.01)
			stand 1 0.02(0.00) 0.14(0.00) -0.05(0.00) -0.17(0.00) 10.27(0.00) 9.66(0.00)
		10	R1 2	0.06(0.00) 0.14(0.00) -0.23(0.00) -0.49(0.00) 12.47(0.01) 11.53(0.01)
			R1 1	0.04(0.00) 0.15(0.00) -0.05(0.00) -0.19(0.00) 10.87(0.01) 9.79(0.00)
			a priori	0(0.00)	0(0.00)	0(0.00)	0(0.00)	10(0.00)	10(0.00)
			stand 2 0.08(0.02) 0.37(0.01) 0.02(0.01) -0.21(0.00) 17.24(0.01) 17.25(0.02)
			stand 1 0.02(0.01) 0.15(0.00) 0.06(0.01) 0.06(0.00) 14.57(0.01) 15.02(0.00)
	with	15	R1 2	0.05(0.02) 0.31(0.01) 0.01(0.01) -0.21(0.00) 18.01(0.02) 17.70(0.01)
			R1 1	0.02(0.01) 0.16(0.00) 0.08(0.00) 0.03(0.00) 15.20(0.01) 15.53(0.00)
			a priori	0(0.00)	0(0.00)	0(0.00)	0(0.00)	15(0.00)	15(0.00)
			stand 2 0.06(0.04) 0.19(0.01) 0.16(0.04) -0.01(0.01) 21.63(0.01) 21.13(0.00)
			stand 1 0.03(0.03) 0.00(0.00) 0.23(0.03) 0.22(0.01) 18.08(0.01) 17.65(0.01)
		20	R1 2	0.06(0.04) 0.16(0.00) 0.14(0.04) -0.03(0.01) 22.61(0.02) 21.86(0.00)
			R1 1	0.02(0.03) 0.01(0.00) 0.25(0.03) 0.20(0.01) 18.55(0.02) 17.94(0.01)
			a priori	0(0.00)	0(0.00)	0(0.00)	0(0.00)	20(0.00)	20(0.00)
					(b) with noise

  (λ 750 nm ) (cm -1 ) 0.12 0.16 0.19 0.23 0.27 0.30 0.34 0.38 0.42 0.45 0.49

	Inclusion	1	2	3	4	5	6	7	8	9	10	11
	Hb (µM * )	25	36	47	58	69	80	91	102 113 125 135
	HbO 2 (µM * )	32	34	36	38	40	42	43	45	47	49	51
	µ a

* µM = µMolar = µmol/L.

Table 2 . 3 :

 23 Summary of the 11 simulated phantoms. equivalent to a composition of 25S λ Hb + 32S λ HbO 2 where S λ Hb and S λ HbO 2 are the reference molar extinction coefficients of Hb and HbO 2 depending on the wavelength λ (Figure 1.1) and an inclusion. The inclusion embedded at 10 mm depth is alternatively one of the 11 different inclusions with gradual concentrations of the two chromophores given in Table 2.3. The values chosen for the series of phantoms are in the range of medical interest. The spectra of the 11 inclusions are displayed in Figure

Table 3 .

 3

1: Characteristics of experimental IRFs used in this manuscript and provided by the hosted lab (Politecnico di Milano). Non gated (NG) and fast gated (FG) SPAD.

Table 3 .

 3 2.

	Inclusion	size 1 size 2 size 3 size 4 size 5
	Diameter (mm)	3.2	4	5	6.8	8.6
	Height (mm)	3.2	4	5	6.8	8.6
	Volume (mm 3 )	25	50	100	250	500
	Equivalent ∆µ expected a	(over 1 cm 3 sphere) (cm 3 ) 0.056 0.087 0.15	0.37	0.94

Table 3 .

 3

2: Summary of the 5 totally absorbing cylinders of various size with diameters, corresponding volumes, and the equivalent absorption perturbations in an area corresponding of a sphere of 1 cm 3 .

Table 3 .

 3

		ρ		Simulation					Experiment
	depth	(mm)	mode	coefficients of the polynomial fit a b c NL slope (∆x=1) (x=1)	a	coefficients of the polynomial fit b NL slope c (∆x=1) (x=1)
	10 mm	15 30	NG G NG G	-6% 64% 1% -6% 64% 1% -5% 56% 1% -4% 53% 3%	-18% -19% -16% -15%	53% 52% 47% 45%	0% 2% -11% 86% 63% 25% 47% 27% 4% -9% 79% 8%	-1% 9% -26% -24%	62% 51% 63% 60%
	15 mm	15 30	NG G NG G	-4% 35% 5% -4% 39% 2% -4% 36% 1% -4% 36% 3%	-22% -23% -20% -21%	28% 30% 29% 28%	-10% 68% -9% 56% -6% 43% 11% -1% 0% -5% 42% 10%	-29% -31% -27% -26%	48% 38% 32% 31%
	20 mm	15 30	NG G NG G	-1% 14% 6% -2% 19% 3% -2% 20% 2% -2% 20% 3%	-12% -23% -21% -22%	12% 15% 16% 16%	-1% -7% -4% -5%	13% 37% 30% 33%	-6% -2% 4% 2%	-12% -39% -29% -30%	11% 23% 21% 23%

3: Coefficients of the polynomial fit of recovered ∆µ vol a vs. ∆µ expected a as defined in section 3.2.2.5 for Non Gated (NG) and Gated (G) SPAD modes.

Table 3 .

 3

						(∆x=1) (x=1)			c	NL (∆x=1) (x=1) slope
	10 mm	15 30	NG G NG G	0% 101% 0% 101% 0% 100% 0% 100%	0% 0% 0% 0%	-1% 0% 0% 0%	100% 115% -139% 123% 100% 59% -25% 71% 100% -9% 152% 2% 100% -9% 149% 2%	-166% -480% -12% -13%	91% 93% 133% 130%
	15 mm	15 30	NG G NG G	1% 105% -3% 0% 100% 0% 0% 100% 0% 0% 101% 0%	-2% 0% 0% 0%	103% 101% 101% 27% 100% -7% 100% -7%	103% 56% 155% 154%	118% 43% -3% -3%	-197% 96% -9% -9%	99% 110% 142% 140%
	20 mm	15 30	NG G NG G	2% 0% 102% -1% 94% 5% 1% 97% 2% 0% 98% 1%	5% -1% 2% 1%	99% 101% 98% 99%	76% 5% 9% 8%	-144% 177% 111% 49% 108% 34% 112% 33%	-106% 8% 17% 13%	9% 121% 126% 127%

4: [Constrained method] Coefficients of the polynomial fit of recovered ∆µ vol a vs. ∆µ expected a as defined in section 3.2.2.5 for Non Gated (NG) and Gated (G) SPAD modes.

Table 4 .

 4 1: Absorbers tested in water-based and ethanol-based liquid phantoms.

	Absorbers	Form	Solubility	Stability	References	Remarks
	China black ink	liquid	water	weeks	China ink, Rotring
	blue dye				Gouro blue, Gubra
	green dye red dye	powder	water (30 g/L) and ethanol	days to weeks	Verde scuro, Gubra Mogano, Gubra	some fluorescence
	brown dye				Noce scuro, Gubra	reduce dynamic range
	Luwsir dye	dark green powder only water (220 g/L)	hours	LUWSIR3 Luminochem unstable with Intralipid
	cyan jet ink black jet ink magenta jet ink	liquid	water and ethanol	days to weeks	ink bulk for HP ink-jet printer (code: HP364)
	blue food dye red food dye	liquid	water and ethanol	NA	blue E133, red E122, Vahiné	spectra Appendix D
	Lunir dye	dark green powder	acetone (388 g/L) and ethanol	NA	LUNIR5 Luminochem	NIR absorbing dye spectrum Appendix D

Acknowledgements

Abbreviations

2D

Two-Dimensional 3D Three-Dimensional reconstruction of even small absorption perturbations (∆µ a = 0.056 cm -1 ) buried down to 20 mm (see Table 3.5), with good accuracy in terms of depth and lateral resolution. By contrast, the absolute quantification of the perturbation is not accurate and its accuracy decrease upon the increase of the inclusion's volume. Thanks to the encouraging results obtained in the simulations, we have studied the suitability of the SiPM for time-domain diffuse optical tomography with experimental measurements.

x 3.5: Recovered value of x, y and depth localization and percentage error for quantification computed for both simulations and experiments for all inclusions (source-detector distance equals to 30 mm).

Experimental measurements

The effect of the depth (ranging from 10 to 30 mm, columns) for inclusion "size 3" in phantom measurements can be seen in Figure 3. 15(b), where are reported the z-y (first row) and x-y (second row) maps. Up to 25 mm depth, the spatial localization (depth and x-y) are as good as it was for simulations (see Table 3.5 for the values). For the 30 mm depth measurement the perturbation is probably too low to properly detect the inclusion, thus resulting in an unsatisfactory localization of the perturbation. For all depths, there is a larger number of artifacts Green's functions of the reference medium or to find a wavelength where biological tissues are "seen" as homogeneous.

Towards clinical

Currently and for future, our laboratory will pursue clinical projects on specific clinical applications with the MW TR DOT. in particular, the laboratory is implicated in a European network whose focus is the application of diffuse optical methods in monitoring brain injury and trauma. An international team is going to work collectively on this subject in research institutions, hospitals and companies. We can see a different of the shape of these IRFs for the early time.

A.2 Conclusion

It seems that the FWHM, even large, has no effect on the contrast at late time i.e. on late photons. The observed time gap at short time doesn't affect the quantification capacity of deep layers and it is correcting in the reconstruction algorithm.

The main parameter of the IRF which impaired the contrast and thus the quantitative information is the down slope with exponential decay constant, α below 5 ns -1 typically in our case.

Appendix B

Complementary results using non-gated and fast-gated SPADs with totally absorbing objects

This appendix links with section 3.2. We present complementary reconstruction results from the simulation and experimental investigations done on the perturbed phantoms with totally absorbing objects using non-gated and fast-gated SPAD.

First, we compare the detection performances of non-gated and fast-gated SPAD in simulation with a short and large source-detector distance (15 mm and 30 mm). Then we display some of the 3D reconstructions used to obtain the quantified plots of section 3.2 and we confirm the quantification performances.

B.1 Detection performances

Supplementary simulations using non-gated SPADs and fast-gated SPADs with a 15 mm and 30 mm source-detector distance were generated to show the interest of gated detectors and to compare with SiPM in section 3.3.3.3.

B.2 Quantification performances

In this section, we show only the 3D reconstructions with the 5 inclusions at 10 mm depth using the fast-gated SPAD in simulation and in experiment because the reconstruction results with the non-gated SPAD were already presented in section 3. For all the z-y slices, we get an accurate localization in depth at 10 mm in z and the red spots almost fit the 1 cm 3 sphere volume represented by a black circle. This confirms the equivalence relation defined in [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF] that a totally absorbing object of 0.025 to 0.5 cm 3 volume gives a specific perturbation absorption over a 1 cm 3 volume in diffuse optics.

For ρ = 15 mm, we observe more spreading problems below the reconstructed inclusion than at ρ = 30 mm. Concerning the µ a quantification in simulation and in experiment, the recovered µ a in the inclusion area increases from inclusion "size 1" to inclusion "size 5" as expected.

In simulation, the retrieved inclusion is more homogeneous than in the experiment where the retrieved µ a shows a peak of absorption. It is probably due to the use of black objects in Concerning the quantification performances using FG SPAD, we are able to distinguish a 0.056 cm -1 absorbing object from a 0.087 cm -1 absorbing object. The increase gradual range of the inclusions is recovered. However the accuracy is difficult to obtain after 0.15 cm -1 and an important apparent decrease of the reconstructed µ a values happens when the inclusion is moved deeper than 10 mm depth. We have demonstrated in the study with SPAD in section 3.2 that with a constrained method adding spatial a priori information, we improve significantly the quantification at all the depths and we achieve accurate reconstructed µ a 3D maps.

Appendix C

Reconstruction results using SiPMs with a 15 mm source-detector distance probe

This Appendix links with section 3.3. It shows the results using the SiPMs with a 15 mm interfiber distance probe in simulation and in experiment. First, we present the detection performances of SiPM at this source-detector distance and then the quantification performances with different inclusions at 10 mm depth. For the inclusions, we use the totally absorbing objects defined in [START_REF] Martelli | Phantoms for diffuse optical imaging based on totally absorbing objects, part 1: basic concepts[END_REF] and explained in section 3.2.2.3. The equivalent absorption perturbations are given in Table 3.2.

C.1 Detection performances

.1 displays the slices of the 3D reconstructed maps of the inclusion " size 3" which has an equivalent ∆µ a (over 1 cm 3 ) equal to 0.15 cm -1 (see Table 3.2). From left to right, the slices show the detection capacity of the inclusion at different depths from 10 to 30 mm with a step of 5 mm. We can distinguish a small red spot representing the perturbation object. Compared to the Figure 3.15 with a large source-detector distance, the detection begins to be difficult for smaller depth. column. We distinguish a red spot surrounding by an almost homogeneous medium of around 0.07 cm -1 in absorption. For ρ = 30 mm, artefacts form a kind of right cross in the x-y directions which is probably due to the scanning pattern. We observe the same thing with the SPAD results (see Appendix B). Some artefacts in the x-y directions are similar in simulation and in experiment which results probably from the scanning pattern and the discretized mesh chosen for the reconstruction process.

C.2 Quantification performances

Concerning the quantification, if we compare with the SPAD results in Appendix B, we identify the same increase of absolute absorption coefficient according to the inclusions. It is confirmed by the quantified plots in Figure 3.18 where we can see the limit of quantification with our reconstruction method.

C.3 Conclusion

This appendix shows the results with a 15 mm interfiber distance probe for SiPM detectors.

The detection performances are slightly reduced compared to the results with the 30 mm sourcedetector distance presented in section 3.3.3. We can well detect up to 20 mm depth and the x and y localization is accurately determined with less than 1 mm of error.

The quantification performances are the same as the large source-detector distance. We get bad results for high absorbing objects and for deep objects ( 15 mm depth). If we consider an error of less than 25% as an accurate quantified result, we can say that the SiPM enables correct quantification up to 0.15 cm -1 at 10 mm depth as with SPAD detector.

D.1 Spectra of food dyes

We diluted red and blue food dyes in water and in ethanol. We measured their absorption spectra with a spectrophotometer (Cary 300 UV-visible, Agilent).

E.1 Protocol for PDMS inclusions

In Chapter 4, in order to assess the quantitation performances and to draw gradual quantified plots of absorbing heterogeneities, we realized 5 gradual bi-absorber inclusion compositions and one inclusion with the same solution as the background phantom. We used different types of phantoms to built the inclusions. In this Appendix, we present the protocol to built the PDMS inclusions.

The objective was to prepare a series of inclusions, all with µ ′ s (600 nm) = 13.5 cm -1 and various µ a with various fractions of black ink and cyan ink (reference: ink bulk for HP ink-jet printer (HP364)). A first inclusion contains only cyan ink and should have µ a (600 nm) = 0.095 cm -1 like the background. Then the series starts with an inclusion containing 100% black ink with and we extrapolate it for our study. The 100% black spectrum is in accordance which what we want. For the other spectra, we have a small offset of around 2 cm -1 but we keep the same curve shape. So, the reduced scattering coefficients of the phantoms are close to the theoretical values. Concerning the absorption coefficients, Figure E.2(a) shows the spectra of each PDMS phantoms represented by the continuous lines. We see that in the part of the visible range (550 to 620 nm) we can slightly distinguish the five spectra associated to the PDMS inclusions and they are in the correct order as a function of their respective fractions of inks. However in the NIR range, the spectra have the same shape. The shapes are significantly different in comparison with what we get in water. The black and cyan curves intersect at around 600 nm with µ a = 0.26 cm -1 which is close to what we wish (µ a = 0.3 cm -1 ). In parallel with the measurements with the spectrophotometer, the shapes of the spectra have the same behaviors but they are flatter with the time-resolved spectroscopy method probably due to the effect of the reduced scattering coefficient and the errors on its estimation. So, the PDMS material flattens the optical signatures of black and cyan jet inks. So, other absorbers have to be tested with this material to be useful for DOT phantoms.

The advantages of this base material are that scattering and absorption coefficients are tunable, the protocol is optimal to obtain the correct µ ′ s and the good proportion of µ a but absorption accuracy is hard to obtain. Furthermore it was demonstrated that the PDMS stiff matrix modifies the absorption spectra. The main manufacturing problems are the estimation of µ a and to have a homogeneous medium without T iO 2 sediments.

Appendix F Preliminary multispectral results on bi-inclusions phantoms using SiPMs

In this Appendix, we present the first results of the experimental campaign on bi-inclusions phantoms at a distance of 25 mm from their centroids (14 mm from their borders). These measurements were carried out very recently in collaboration with the Politecnico di Milano. With the same protocol than the experiments with the "balloons" inclusions in Chapter 4, we put in the same liquid phantom: two "balloons" inclusions of same size (1 mL) and same composition (50% cyan and 50% black inks), two "balloons" inclusions of different size (1 mL and 2 mL) but of same composition (50% cyan and 50% black inks) and two "balloons" inclusions of same size (1 mL) and different composition (100% cyan ink and 50% cyan and 50% black inks) all embedded at different depths (see Picture F.1(a)). We moved alternatively each inclusion The data processing was the same as previously but we increased the size of the scanning mesh in the x-direction. Figure F.2 shows preliminary 3D reconstructions with the 550-580-600-nm wavelength combination of a pair of balloons of the same size and different composition ("100% cyan" and "50% cyan 50% black") that we placed one after the other from 10 to 15 mm depth.

The first row shows the z-x slices for y = 0 where can be seen a yellow spot for the left balloon on cyan representations and inversely a yellow spot for the right balloon on black representations for the two first columns. In the second row with z = 10 mm and the last column, the right balloon has more quantity of cyan ink than the left balloon and the inversion is obtained at z = 15 mm because the 100% cyan inclusion (left) is embedded at 15 mm depth. We can well identify and localize the two inclusions in these three cases. In the last row with z = 20 mm, we reconstruct only the background with some fluctuations (+/-10%) with around 33% for cyan and 0% for black inks as expected.