1. Introduction

Near-infrared optical measurements in biological tissues offer relevant contrast on phenomena related to hemoglobin and other endogenous chromophores like water or lipids. Tumors may be detected because of a tissue composition differing from surrounding healthy tissue as in mammography [1], brain activity can be observed due to changes in oxy and deoxyhemoglobin concentrations [2,3], and hemorrhages can be detected in the brain of a premature infant [4]. For some applications, detection is the final goal, but in other cases like biopsy guidance, precise localization and quantification of the contrast are the clinically relevant pieces of information.

In most cases for human applications, transmission measurements cannot be obtained, so only reflectance measurements are available to retrieve information on optical contrast in depth. For the particular case of prostate measurements, the endorectal probe has limited distances between sources and detectors [5]. Even in cases where transmission measurements could be possible due to anatomy, probe-like configurations enabling only reflectance measurements with small source-detector distances can be preferred. Indeed, reflectance-only devices are more portable, thus facilitating medical practice in applications such as mammography [6] and offer a better control on distances between sources and detectors which is interesting for brain applications, where the use of remote optical fibers in flexible helmets is a cause of localization errors [7].

We therefore investigate the probe-like geometry for reflectance measurements characterized by small distances between sources and detectors. The particular challenge of this configuration is to retrieve information of deep layers. This information is known to be carried by the so-called late photons, which are much less prevalent than first photons coming back from superficial layers [8]. So, continuous-wave or frequency domain setups are not adapted to image deep layers since the very few late photons are shadowed by the much larger number of first photons. Time-resolved acquisition chains are intrinsically more appropriate to detect late photons with high signal to noise ratio.

In the past years, multiple solutions have been proposed and studied for processing time-resolved optical signals in order to maximize the use of their information content, for fluorescence [9–11] and endogenous diffuse optical imaging [12–14]. In a previous publication [15], we have exposed an analysis method based on the Mellin-Laplace transform to process time-resolved signals for Diffuse Optical Tomography (DOT). We have shown that this approach enables time-windowing of the diffusion curve and detection of deep absorbing inclusions. This method allows to tune the sharpness of the analysis of diffusion curves, thus provides an intermediate scheme between the moments [13] and full time-resolved analysis [14].

In the present article, we aim at showing more precisely how the Mellin-Laplace transform (MLT) can be applied to process experimentally acquired time-resolved signals, so called Time-Point Spread Functions (TPSFs). For that purpose, we have developed an original way to take into account the Instrument Response Function (IRF) by using a reference measurement in a known diffusive medium. We also show how the choice of MLT orders considered in the analysis is important for a proper detection and localization of an absorbing inclusion. From this, we investigate the performance of the method in terms of detection and localization of deep absorbing inclusions with reflectance measurements, depending on the dynamic range of the acquired or simulated TPSFs, referred to as “signal dynamics”. In this article, the method is first illustrated on simulated signals featuring the most important sources of experimental noise and then on real measurements.

This article is organized as follows. In the first place, the reconstruction method is developed. Only the aspects dealing specifically with experimental signals are developed here, the other aspects of the method being detailed in [15]. In a second chapter, the method is illustrated on simulated signals for the “small” interfiber distance of 1.5 cm. Optimization with respect to signal dynamics and performance in terms of detection and localization are discussed. To finish, the applicability of this method to experimentally acquired signals is demonstrated on measurements obtained with a time-resolved acquisition setup including a femtosecond laser, photomultipliers and a time-correlated single photon counting (TCSPC) board.

2. Reconstruction method

2.1. DOT with the Mellin-Laplace transform

The followed approach for processing time-domain DOT with the Mellin-Laplace transform (MLT) is fully detailed in [15]. The method was developed considering the formalism of the diffusion approximation to describe light propagation in a turbid medium in which the absorption coefficient is small with respect to the reduced scattering coefficient (µ_a << µ_s′) [16]. In this context, the perturbation of time-domain measurements due to small variations of the absorption coefficient in the medium is expressed as follows (* is the convolution operator with respect to the time variable):

where $G_i^{\beta }(\overrightarrow{r},t)$ and $G_i^{\alpha }(\overrightarrow{r},t)$ are the “Green's functions” solutions of the time-resolved diffusion equation with the source terms $\delta (\overrightarrow{r}-{\overrightarrow{r}}_i)\delta (t)$, i = s or d relates to the position of the source and detector respectively, and α or β are two different media. These media differ by the absorption coefficient:$\delta {\mu }_a(\overrightarrow{r})={\mu }_a^{\beta }(\overrightarrow{r})-{\mu }_a^{\alpha }(\overrightarrow{r})$.

Equation (1) is the fundamental relation on which iterative DOT reconstruction algorithms are based [17].

The Mellin-Laplace transform M^(p,n) as defined in Eq. (2) can be applied to a convolution product of multiple time-resolved functions and calculated without explicit calculation of this convolution product (Eq. (3)):

where p (in ns⁻¹) is a real positive and n is a positive or null integer,

The property of Eq. (3) enables the computation of Eq. (1) without discretization into small time steps dt and without explicit calculation of the convolution product. The analysis can be carried out with a tunable precision: 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. In [15], reconstructions of µ_a maps using this scheme were obtained considering an ideal IRF equal to a Dirac function$\delta (t)$.

2.2. Treatment of the IRF

The perturbation approach described in Eq. (1) links $\delta {\mu }_a(\overrightarrow{r})$ to variations between the Green's functions of two media. However, the Green's functions are not accessible experimentally: the measurements correspond to these functions convolved by the IRF. In order to approach $\delta {\mu }_a(\overrightarrow{r})$ based on variations between measurements acquired in two media, one needs to take into account the real IRF. Indeed, the laser sources used for DOT have a time spread (from a few hundreds of femtoseconds to a few picoseconds), the detectors like photomultipliers tubes (PMT) or single-photon avalanche diodes (SPAD) have typical responses in the order of a few hundreds picoseconds full width at half maximum, and the optical fibers guiding the light cause time-spread of the signal as well. Moreover, these behaviors can be slightly different depending on each detector and fiber and change in time with instrumental drifts (mainly due to temperature). The real IRF needs to be accounted for in the reconstruction scheme: it needs to be considered separately for each couple of source and detector and calibrated as often as possible to account for the time drifts.

State of the art. This issue has been largely studied in the past ten years and many strategies were developed. One solution is to measure directly the IRF by placing each source in front of each detector but it has the drawback of being time-consuming [18]. A second solution is to calibrate in situ the IRF by measuring its reflection on the surface of the object under study [19]. This solution was performed with success in vivo on the head of newborn infants [3]. Another solution is to use a reference phantom to calibrate the IRF, so as to work with relative and not absolute datatypes [20]. This solution has two requirements: the reference phantom should have the same geometry as the object to study but it should also have the same optical properties as its background so that the use of relative datatypes is allowed.

Proposed approach. Our developed approach requires a measurement on a reference phantom with known properties, but these properties do not have to be exactly similar to the background of the object under study, as the analysis is not based on relative datatypes. The novelty of the method consists in applying the MLT on a combination of the reference measurement and the measurement of the object, without having to explicitly deconvolve the IRF at any time.

The reconstruction process requires the TPSFs $A_{sd}^{}(t)$measured for all couples of sources s and detectors d in a reference medium A and the TPSFs $B_{sd}^{}(t)$measured in the medium B under study. The whole iterative process is summarized in Fig. 1 .

 

Fig. 1 Schematic view of the followed algorithm for tomographic reconstruction.

Download Full Size | PDF

The crucial relations and steps of this process are the following. If $G_{sd}^A(t)$and $G_{sd}^B(t)$ are the ideal TPSFs in media A and B for Dirac sources and for couples of sources s and detectors d, the measured TPSFs are their time convolution with the IRF for each couple of source s and detector d:

Due to convolution properties, the following equation can then be deduced from Eq. (4):

An iterative process can be carried out to reconstruct the absorption coefficients ${\mu }_a^B(\overrightarrow{r})$of all points r in medium B (See Fig. 1) so that $G_{sd}^B(t)$verifies Eq. (5). Let ${\mu }_{a,k}^B(\overrightarrow{r})$be the absorption coefficient map in medium B estimated at iteration step k. With this map we solve numerically the diffusion equation to get the Green functions$G_{sd,k}^B(t)$. Then, the true TPSFs with the real set of coefficients ${\mu }_a^B(\overrightarrow{r})$ can be approached by the estimation$G_{sd,k}^B(t)$ described in Eq. (1):

with $\delta {\mu }_{a,k}(\overrightarrow{r})={\mu }_a^B(\overrightarrow{r})-{\mu }_{a,k}^B(\overrightarrow{r})$ and $G_{s,k}^B(\overrightarrow{r},t)$ and $G_{d,k}^B(\overrightarrow{r},t)$ respectively the theoretical TPSFs from source s or detector d to a point r in medium B, at iteration k. At iteration k, we can therefore link the measures and the estimated TPSFs $G_{sd,k}^B(t)$ with variations in the absorption map $\delta {\mu }_{a,k}(\overrightarrow{r})$ by combining Eq. (5) and Eq. (6):

and get ${\mu }_{a,k+1}(\overrightarrow{r})={\mu }_{a,k}(\overrightarrow{r})+\delta {\mu }_{a,k}(\overrightarrow{r})$. µ_a is therefore nonlinearly reconstructed by this iterative algorithm.

Instead of discretizing Eq. (7) with time steps dt which would be computationally intense, we propose to calculate its MLT (Eq. (2)). By applying the time convolution properties of the MLT (Eq. (3)) to Eq. (7), we obtain the following linear system:

Equation (8) has the form of a linear matrix equation Y = WX where X contains the$\delta {\mu }_{a,k}(\overrightarrow{r})$that yields an estimation of the map of absorption coefficients ${\mu }_a^B(\overrightarrow{r})$ in medium B. At the iteration step k, the MLT of the measurements in the reference medium $A_{sd}^{(p,e)}$ and in the medium to study $B_{sd}^{(p,e)}$are therefore linked to the MLT of the computed Green's functions $G_{sd,k}^{A(p,f)}$and$G_{sd,k}^{B(p,f)}$ in these media at this iteration, and to$\delta {\mu }_{a,k}(\overrightarrow{r})$. Vector Y is constructed by concatenating all MLT orders (N_MLT in total) for each couple of source (N_s in total) and detector (N_d in total), its dimension is $(N_{MLT}\times N_s\times N_d)\times 1$. The Jacobian matrix W is obtained similarly and has a dimension of $(N_{MLT}\times N_s\times N_d)\times N_x$ with N_x the number of pixels defined in the medium.

2.3. Solving the inverse problem

The inverse problem is solved as described in [15]. At each step k we minimize ${\Vert L\cdot W\cdot R\cdot X^{\prime }-L\cdot Y\Vert }^2$using a gradient descent. R is the right preconditioner. We chose a diagonal matrix whose elements R_mm correspond to the square of the distance between the pixel m and the set of source and detector locations. The left preconditioner L is a diagonal matrix whose diagonal elements estimate the variance on each order of the MLT of Y. We then obtain X by a matrix multiplication by R: $X=R\cdot X^{\prime }$.

We notice that convergence is reached after less than 10 iterations in all the tested cases (the criteria is that the reconstructed depth of the absorbing inclusion is stable). All the results presented in this article are taken from the 15th iteration.

2.4. Analysis of the reconstructed images

Objective criteria are calculated from the reconstructed images in order to assess the performance of the method. These criteria include localization of the inclusion in x and z directions and reconstructed µ_a mean value.

From the reconstructed image (Fig. 2(a) ), we consider as belonging to the absorbing inclusion the pixels whose δµ_a value with respect to the background is at least 50% of the maximum value of δµ_a in the image (Fig. 2(b)). We chose the threshold of 50% because the reconstructed absorbing spots are very large, due to the choice of minimizing the 2-norm. We then calculate the x and z position of the center of mass of the absorption coefficient in this 50% spot and the mean value of the reconstructed µ_a in this spot (Fig. 2(b)).

 

Fig. 2 (a) An example of 2D reconstructed map of µ_a in the medium. (b) The 50% spot and the extracted parameters: x and z position of its center of mass and the mean µ_a in this spot.

Download Full Size | PDF

3. Simulations

3.1. Numerical phantom and probe

Phantom. We present here results obtained with a code in two dimensions. To simulate the studied medium, we use a two-dimensional numerical phantom of 6x5 cm, with a rectangular mesh grid and mesh steps of 0.1 cm in x and z directions (Fig. 3 ). We implement the extrapolated boundary conditions as follows: we enlarge the medium by 1 mm and impose a null photon density at the extrapolated boundary. The absorption inclusion is a disc (diameter set to 1 cm for the whole study) placed at the center of the medium in the x direction. Depth of this inclusion, defined as the distance between the surface of the medium and the center of this disc, varies from 10 to 35 mm with steps of 5 mm. The reduced scattering coefficient of the background and the inclusion are set to μ′_s = 10 cm⁻¹ for the whole study, a typical value for biological tissues in the near-infrared window [21]. The background of the medium is simulating the breast with µ_a = 0.1 cm⁻¹. The inclusion has the following optical properties: µ_a = 0.6 cm⁻¹ and μ′_s = 10 cm⁻¹.

 

Fig. 3 Numerical 2D phantom and probe (crosses: sources “S”, circles: detectors “D”, large grey disc: absorption inclusion). The mesh is represented only in 1 cm². Only the couples S1D2 and S2D1 are used for the analysis.

Download Full Size | PDF

Probe. The numerical probe is similar to the experimental one (Subsection 4.1) and made of 2 excitation fibers and 2 detection fibers. The (x, z) coordinates in mm of the source points are [S1 (−5, −1); S2 (+5, −1)] and of the detection points are [D1 (−10, 0); D2 (+10, 0)]. To comply with the hypothesis of the isotropic point source, the sources are placed at 1/μ′_s = 0.1 cm inside the medium (Fig. 3).

For the following analysis, only the couples S1D2 and S2D1, distant by 1.5 cm will be used. These two couples will enable the localization of the inclusion in the x direction below the probe. The inclusion being placed exactly in the middle of the two sources, the contrast measured by the couples S1D2 and S2D1 should be similar unless noise becomes preponderant, which should cause an error in the x localization of the inclusion. With this compact geometry, the probe has a width of 2 cm only, includes interfiber distances of 1.5 cm and enables localization in the x direction.

3.2. Simulation of signals

Generation of TPSFs. The reference measurement is taken in the homogeneous medium without the inclusion and the other measurement is taken in the same medium with the inclusion. Signals are simulated by solving the time-resolved diffusion equation for Dirac point sources with the finite-volume method in the 2D grid of the numerical phantom (implemented in Matlab®) for the medium with and without the inclusion. One TPSF is simulated per couple of source and detector and per depth of the inclusion (6 values of depth in total). To insure good accuracy of the computation, the time step is dt = 3 ps, which also corresponds to the time resolution of our TCSPC setup (Subsection 4.1). Thus, our simulated signals have the same time sampling as experimental ones.

Adding noise. Real experimental signals are degraded by two main sources of noise: the Dark Count Rate (DCR) and the photonic noise or shot noise. The DCR is due to thermal photogeneration of charges in the detector: on the measured TPSFs, it concretely consists in an offset present in all time channels. The value of this offset depends on the type of detector and temperature. The term “photonic noise” describes the fact that detectable photoelectrons show a statistical variation. Photonic noise is modeled by a Poisson distribution: its variance is equal to the number of counted photons, so this source of noise becomes less predominant when this number increases, i.e. when acquisition time increases.

To add these two sources of noise to our simulated TPSFs, we first convert them to a number of photons per channel by multiplying them by a constant. This constant is calculated for each TPSF to obtain the desired total number of photons which will define the signal dynamics. Here, we will consider alternatively 10⁶ and 10⁸ photons per TPSF, which correspond to realistic dynamics encountered on experimental acquisitions. An offset of 100 photons per time channels is added to simulate DCR. The last step is to add Poisson noise to each time channel. We repeat this last step 10 times to simulate a set of 10 measurements for each TPSF (the reference and each depth of the inclusion in the heterogeneous case), in order to evaluate the sensitivity of our technique to photonic noise.

DCR correction. Before processing the simulated signals, we remove the offset caused by DCR, calculated on a flat portion of the signal (between 6 and 8 ns). The corrected signals are shown in Fig. 4(a) and Fig. 4(b). Only time channels belonging to the diffusion curve are kept in the analysis: channels from 0 to 2.7 ns for 10⁶ photons and 0 to 4.3 ns for 10⁸ photons.

 

Fig. 4 Simulated signals without the inclusion (red) and with the inclusion at 2 cm depth (blue), (a) 10⁶ photons after DCR correction, (b) 10⁸ photons after DCR correction. For all: 1 noise realization only.

Download Full Size | PDF

Variance. We evaluate the variance on each time channel as the number of photons per time-channel before the DCR correction.

3.3. Results

Contrast analysis. As we considered as reference the medium to study without the inclusion, we use the MLT orders in the reference and in the studied medium to calculate a contrast (%):

The analysis and reconstruction results presented in this article were carried out for a given precision of MLT of p = 3 ns⁻¹, which represents a good compromise between the pieces of information extracted from the TPSF and the computation time.

Figure 5 depicts the contrast on the MLT orders calculated with Eq. (2) for p = 3 ns⁻¹, for each depth of the inclusion and one couple of source and detector, both for 10⁶ and 10⁸ photons. The results with the other couple of source and detector are not displayed as they are identical, the inclusion being symmetric with respect to these two couples. For all values of depth, the contrast is low for the first orders related to first photons and increases with the following orders related to later photons. Error bars increase as well when the order increases. Additionally, improving the signal dynamics enables to decrease the error bars: for a given order of MLT, the error bar decreases when going from 10⁶ to 10⁸ photons.

 

Fig. 5 Contrast on MLT orders (p = 3 ns⁻¹) for the different depths for S1D2, (a) 10⁶ photons, (b) 10⁸ photons (error bars are obtained from the 10 noise realizations).

Download Full Size | PDF

Considering now each depth separately on Fig. 5, we clearly see that the deeper the inclusion, the lower the contrast. For 10⁶ photons, inclusions at depths between 10 and 25 mm give contrast values clearly above the noise level, for orders 0 to 30 for 10, 15 and 20 mm and only for orders 0 to 20 for 25 mm. At depths larger than 30 mm, the contrast due to the inclusion is not statistically relevant, so we cannot expect to detect these inclusions. For 10⁸ photons, the contrast is above the noise level up to 35 mm included. Figure 5 therefore suggests that the contrast is statistically relevant on higher orders of MLT when using a higher signal dynamics.

Reconstruction results. The 2D maps of µ_a are reconstructed on the same grid as used for the direct model and depicted in Fig. 3. As it is not possible to present all reconstructed images, the obtained results are summarized in Fig. 6 . For each signal dynamics and each depth of the inclusion, separated reconstructions were carried out for the 10 noise realizations. We analyzed the resulting images with the criteria described in Subsection 2.4. We tested different choices of maximum MLT orders for solving the inverse problem. The importance of this maximum order is revealed. For a given precision p, the quality of the reconstructions is very sensitive to the choice of the maximum orders of MLT included for solving the inverse problem. Whereas the optimum choice enables the best detection and localization of the inclusion allowed by the algorithm, including too few orders causes low sensitivity to deep inclusions and strong underestimation of depth. High orders are more subject to photonic noise so that their inclusion in the reconstruction degrades the results.

 

Fig. 6 Reconstructed depth versus true depth for different choices of maximum orders of MLT included for the reconstruction, (a) 10⁶ photons, (b) 10⁸ photons (error bars obtained from the 10 noise realizations).

Download Full Size | PDF

For example, as it can be seen in Fig. 6, for 10⁶ photons, including only orders n = 1 to 5 in the analysis causes a large underestimation of depth from 10 to 25 mm. Adding orders up to n = 15 yields a more accurate estimation of the depth. However, adding even more orders (up to n = 30) degrades again the localization of the inclusion (underestimation and increased error bars). The phenomenon is similar with 10⁸ photons, were the optimum maximum order to include is n = 30. These results are consistent with the analysis of the contrast graph of Fig. 5.

Reconstructed images for one noise realization are presented in Fig. 7 . They were obtained for the optimum configurations deduced from Fig. 6, including orders n = 1 to 15 for 10⁶ photons and n = 1 to 30 for 10⁸ photons. These images enable to see the limit of the proposed method depending on the signal dynamics that is given by the number of detected photons. In the range 10 mm to 15 mm, reconstructed images are similar for both signal dynamics. From 20 mm to 25 mm, the inclusion is detected for both cases but the estimated value of µ_a and depth are more accurate with 10⁸ photons. From 30 mm to 35 mm, only the larger dynamics enables to detect the inclusion but its depth is underestimated, as we could already see in Fig. 6. Let us note as well that the deeper the inclusion, the larger the size of the reconstructed spot and the lower the estimated µ_a. This underestimation of µ_a in depth could be limited by improving the spatial resolution of the image, for example by multiplying the number of sources and detectors or at the algorithm level by using priors or sparsity.

 

Fig. 7 Reconstructed images of µ_a for 10⁶ photons and 10⁸ photons (one noise realization). The red dotted circle depicts the true position and size of the absorbing inclusion.

Download Full Size | PDF

Figure 7 depicts the results obtained for one noise realization so it cannot reflect the variability of results due to photonic noise. Figure 8 summarizes the performance of the method for the two studied signal dynamics, with error bars obtained from the 10 noise realizations. For all studied parameters (depth, x and µ_a) and both signal dynamics, error bars increase with depth. They are however always larger for 10⁶ photons than for 10⁸ photons. For 10⁶ photons, the mean reconstructed depth does not increase after 25 mm: reconstructions are driven by photonic noise and not by the presence of the inclusion. From 15 mm to deeper positions, the depth is always better estimated with 10⁸ photons than with 10⁶ photons. Increasing the signal dynamics enables not only to detect deeper inclusions but also to estimate better their depth. However, from 20 mm to deeper, the depth is underestimated even with 10⁸ photons. This effect was already observed by Selb et al. starting from the depth of 25 mm for an interfiber distance of 25 mm [22].

 

Fig. 8 Summarized reconstruction performances for 10⁶ and 10⁸ photons for 10 noise realizations, (a) reconstructed depth, (b) x localization, (c) Mean reconstructed µ_a.

Download Full Size | PDF

4. Experiments

4.1. Setup and protocol

The measurement setup put in place and depicted in Fig. 9 is similar to time-resolved imaging systems suggested by multiple sources in literature, like Kacprzak et al. [23].The excitation source consists of a Ti:sapphire femtosecond laser (Chameleon, Coherent, Santa Clara, USA) operated at 80 MHz and 770 nm injected in an optical fiber. A mean power of 5 mW is delivered by each excitation fiber. Light is collected with two detection fibers, each of them being directed to a photomultiplier tube (PMT) (PMC-100-20, Becker&Hickl, Germany) connected to one channel of a TCSPC board (SPC-134, Becker&Hickl, Germany). The IRF measured in this setup is 190 ps FWHM. Optical attenuators are placed in front of the PMT in order to reach exactly the limit count rate of 10⁶ photons per second.

 

Fig. 9 TCSPC setup, probe and optical phantom.

Download Full Size | PDF

The probe has the same geometry as described in Subsection 3.1. The background medium of the optical phantom is made of Intralipid^© and black ink (Rotring, Art. 591 017) to reach the same properties as the simulated phantom (µ_a = 0.1 cm⁻¹ and µ_s′ = 10 cm⁻¹). A solid absorbing inclusion consisting of a cylinder of 8 mm in diameter and 6 cm long made of resin, TiO₂ particles and the same black ink, is placed in the background medium at the center of S1 and S2. Its optical properties are µ_a = 0.6 cm⁻¹ and µ′_s = 10 cm⁻¹. During the experiment, the inclusion is moved at different depths: from 10 mm to 35 mm from the surface with steps of 5 mm. Before and after each of these “depth” measurements a “reference” is acquired in the medium without the inclusion. The displacement of the inclusion in the z direction is automated with a translation stage. For each measurement (“depth” or “reference”), acquisitions of 1 second are performed successively for the two couples of source and detector.

4.2. Results

Figure 10 shows all the “depth” and “reference” raw time-resolved signals acquired for S1D2. Figure 10(a) shows that all the measured references are superimposed, which assesses that no major drift was observed during the whole experiment. On the raw “depth” signals of Fig. 10(b), we can easily distinguish the measurements at the depths of 10 mm and 15 mm from the others, which seem to be superimposed.

 

Fig. 10 Raw signals (after DCR correction) for S1D2, (a) All references, (b) All “depth” measurements (depth 1 acquired between refs [1] and [2], etc.).

Download Full Size | PDF

Figure 11 depicts the contrasts calculated with Eq. (9) for different orders of MLT with p = 3 ns⁻¹, on both couples of source and detector. The error bars were calculated from the contrasts obtained when using the reference before and after the “depth” measurement. They indicate the sensitivity of the contrast to photonic noise. The values of contrast for the depths of 10 mm and 15 mm are very high. The depth of 20 mm shows a visible contrast as well. However, from 25 mm to 35 mm the contrast is null. We can notice on Fig. 11 that the contrast profiles obtained with experimental signals are different from those obtained on simulations with an ideal Dirac IRF (Fig. 5). This is due to the fact that the IRF of our setup features a “tail”, which tends to broaden the signals. Other types of detectors (like single-photon avalanche diodes or hybrid PMTs), optimization of the optical path (optical fibers, etc.) could be considered to obtain a thinner IRF. The impact of different non ideal IRF on the reconstructions with our method is still to be studied.

 

Fig. 11 Contrast on MLT orders for different depths of the absorbing inclusion.

Download Full Size | PDF

For each depth value, we ran two reconstructions on the same mesh grid as described in Fig. 3. One reconstruction was realized using the reference acquired before the “depth” measurement and another one using the reference acquired after. This was done in order to evaluate the sensitivity of the reconstruction to photonic noise. Figure 12 synthesizes the obtained results (MLT parameters: p = 3 ns⁻¹, orders n = 1 to 17).

 

Fig. 12 Summarized reconstruction performances for experimental signals with an acquisition time of 1 s, (a) reconstructed depth, (b) x localization, (c) Mean reconstructed µ_a.

Download Full Size | PDF

Figure 13 shows that from 10 mm to 20 mm, the inclusion is well detected and reconstructed at the proper depth. The results are not influenced by photonic noise, as the images obtained with two different reference measurements are similar. From 25 mm to 35 mm, the obtained images seem to be only resulting from photonic noise: the reconstructed µ_a is very low for all of them and the results are very different depending on the chosen reference. These results are consistent with the order of magnitudes obtained on simulations: with 10⁶ photons we could see that the depth limit in detection was at 25 mm for a similar contrast of the inclusion.

 

Fig. 13 Reconstructed images at iteration 15 for all depths of the inclusion, for 2 different references. From 10 to 20 mm, the reconstructed images are similar for both references. From 25 mm to deeper, analogous reconstructions obtained with the 2 different references differ very much: they are driven by photonic noise.

Download Full Size | PDF

5. Discussion

The results obtained in simulations and experiments indicate that the proposed method enables to extract the information of late photons to detect and localize deep absorbing inclusions. This information is included in the MLT of the experimental TPSF for orders n smaller than a given value, depending on the signal dynamics. Detection of deep inclusions is possible only if the signal dynamics is sufficient. However, as mentioned when presenting the experimental setup, acquisition chains like TCSPC are limited with a maximum count rate, which requires long acquisitions in order to increase signal dynamics. It is important to keep in mind that for small interfiber distances the detection and photon counting electronics are responsible of this limitation and not the available source power at the collection fiber. Indeed, the maximum count rate of 10⁶ photons per second is reached for source power well below the safety limits for small interfiber distances [24]. Simulations have shown that in our typical case multiplying by 100 the number of photons (and therefore the acquisition time) enabled to extend the maximum depth detection from 25 mm to 35 mm. This improvement in the imaged depth range is a crucial step to make probes more useful as medical imaging tools. However, such acquisition durations can be prohibitive when transferring this technique to the clinic. Other acquisition modes should be considered to acquire high dynamic signals in minimum time. A promising option is the possibility to combine TCSPC with fast-time gating of the detector to increase the count rate of late photons as it can be done with fast-gated SPAD [24].

In this article, the analysis was carried out for a given precision of the MLT (p = 3 ns⁻¹). Similar results to those exposed here have also been obtained with a higher or smaller precision of the MLT (results not shown here). As stated before, the important point is to restrict the range of MLT orders to those that are significant depending on the signal dynamics. However, literature suggests that a high analysis precision of the TPSF can be interesting for imaging complex objects: Gao et al. have shown improvement of image quality by using full-time resolved data in the cylindrical geometry [14]. Therefore, studies on more complex objects should be carried out to conclude the interest of using high precision.

6. Conclusion

We have shown how to apply the Mellin-Laplace transform as defined in [15] to experimental signals and detailed how to optimize the data analysis procedure. This approach was applied to simulated signals modeled with the main sources of experimental noise. A study on two different signal dynamics has shown how the method should be applied in order to extract the information from the latest photons to detect and localize deep absorbing inclusions. To finish, we have demonstrated the application of this methodology on experimental time-resolved signals. This study shows that this type of probe geometry could be useful to image a biological tissue like the breast in a depth range of 0 to 20 mm or 0 to 35 mm depending on acquisition time. This last point is crucial and currently limited by the detection and not by the available source power. New solutions, like fast-gated SPADs associated to TCSPC could bring a solution to this issue.

Next research steps to follow are to investigate techniques to acquire high signal dynamics in limited acquisition time and to study more complex media.