## Abstract

We show how to apply the Mellin-Laplace transform to process time-resolved reflectance measurements for diffuse optical tomography. We illustrate this method on simulated signals incorporating the main sources of experimental noise and suggest how to fine-tune the method in order to detect the deepest absorbing inclusions and optimize their localization in depth, depending on the dynamic range of the measurement. To finish, we apply this method to measurements acquired with a setup including a femtosecond laser, photomultipliers and a time-correlated single photon counting board. Simulations and experiments are illustrated for a probe featuring the interfiber distance of 1.5 cm and show the potential of time-resolved techniques for imaging absorption contrast in depth with this geometry.

© 2013 OSA

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

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

## 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^{−1} 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^{−1}. The inclusion has the following optical properties: µ_{a} = 0.6 cm^{−1} and μ′_{s} = 10 cm^{−1}.

**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^{6} and 10^{8} 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^{6} photons and 0 to 4.3 ns for 10^{8} photons.

**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^{−1}, 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^{−1}, for each depth of the inclusion and one couple of source and detector, both for 10^{6} and 10^{8} 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^{6} to 10^{8} photons.

Considering now each depth separately on Fig. 5, we clearly see that the deeper the inclusion, the lower the contrast. For 10^{6} 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^{8} 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.

For example, as it can be seen in Fig. 6, for 10^{6} 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^{8} 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^{6} photons and n = 1 to 30 for 10^{8} 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^{8} 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.

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^{6} photons than for 10^{8} photons. For 10^{6} 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^{8} photons than with 10^{6} 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^{8} 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].

## 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^{6} photons per second.

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^{−1} and µ_{s}′ = 10 cm^{−1}). A solid absorbing inclusion consisting of a cylinder of 8 mm in diameter and 6 cm long made of resin, TiO_{2} 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^{−1} and µ′_{s} = 10 cm^{−1}. 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.

Figure 11
depicts the contrasts calculated with Eq. (9) for different orders of MLT with p = 3 ns^{−1}, 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.

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^{−1}, orders n = 1 to 17).

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^{6} photons we could see that the depth limit in detection was at 25 mm for a similar contrast of the inclusion.

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

## References and links

**1. **B. J. Tromberg, B. W. Pogue, K. D. Paulsen, A. G. Yodh, D. A. Boas, and A. E. Cerussi, “Assessing the future of diffuse optical imaging technologies for breast cancer management,” Med. Phys. **35**(6), 2443–2451 (2008). [CrossRef]

**2. **B. Montcel, R. Chabrier, and P. Poulet, “Time-resolved absorption and hemoglobin concentration difference maps: a method to retrieve depth-related information on cerebral hemodynamics,” Opt. Express **14**(25), 12271–12287 (2006). [CrossRef]

**3. **J. C. Hebden, M. Varela, S. Magazov, N. Everdell, A. Gibson, J. Meek, and T. Austin, “Diffuse optical imaging of the newborn infant brain,” in *2012 9th IEEE International Symposium on Biomedical Imaging (ISBI)* (IEEE, 2012), pp. 503–505.

**4. **T. Austin, A. P. Gibson, G. Branco, R. M. Yusof, S. R. Arridge, J. H. Meek, J. S. Wyatt, D. T. Delpy, and J. C. Hebden, “Three dimensional optical imaging of blood volume and oxygenation in the neonatal brain,” Neuroimage **31**(4), 1426–1433 (2006). [CrossRef]

**5. **J. Boutet, L. Herve, M. Debourdeau, L. Guyon, P. Peltie, J.-M. Dinten, L. Saroul, F. Duboeuf, and D. Vray, “Bimodal ultrasound and fluorescence approach for prostate cancer diagnosis,” J. Biomed. Opt. **14**(6), 064001 (2009). [CrossRef]

**6. **J. Liu, A. Li, A. E. Cerussi, and B. J. Tromberg, “Parametric diffuse optical imaging in reflectance geometry,” IEEE J. Quantum Electron. **16**(3), 555–564 (2010). [CrossRef]

**7. **J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, “Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol. **47**(23), 4155–4166 (2002). [CrossRef]

**8. **S. Del Bianco, F. Martelli, and G. Zaccanti, “Penetration depth of light re-emitted by a diffusive medium: theoretical and experimental investigation,” Phys. Med. Biol. **47**(23), 4131–4144 (2002). [CrossRef]

**9. **N. Ducros, L. Hervé, A. Da Silva, J.-M. Dinten, and F. Peyrin, “A comprehensive study of the use of temporal moments in time-resolved diffuse optical tomography: part I. Theoretical material,” Phys. Med. Biol. **54**(23), 7089–7105 (2009). [CrossRef]

**10. **N. Ducros, A. Da Silva, L. Hervé, J.-M. Dinten, and F. Peyrin, “A comprehensive study of the use of temporal moments in time-resolved diffuse optical tomography: part II. Three-dimensional reconstructions,” Phys. Med. Biol. **54**(23), 7107–7119 (2009). [CrossRef]

**11. **J. Riley, M. Hassan, V. Chernomordik, and A. Gandjbakhche, “Choice of data types in time resolved fluorescence enhanced diffuse optical tomography,” Med. Phys. **34**(12), 4890–4900 (2007). [CrossRef]

**12. **H. Zhao, F. Gao, Y. Tanikawa, and Y. Yamada, “Time-resolved diffuse optical tomography and its application to *in vitro* and *in vivo* imaging,” J. Biomed. Opt. **12**(6), 062107 (2007). [CrossRef]

**13. **M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. **44**(7), 1699–1717 (1999). [CrossRef]

**14. **F. Gao, H. Zhao, and Y. Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt. **41**(4), 778–791 (2002). [CrossRef]

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

**16. **A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. **43**(5), 1285–1302 (1998). [CrossRef]

**17. **S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt. **34**(31), 7395–7409 (1995). [CrossRef]

**18. **E. M. C. Hillman, J. C. Hebden, F. E. W. Schmidt, S. R. Arridge, M. Schweiger, H. Dehghani, and D. T. Delpy, “Calibration techniques and datatype extraction for time-resolved optical tomography,” Rev. Sci. Instrum. **71**(9), 3415–3427 (2000). [CrossRef]

**19. **J. C. Hebden, F. M. Gonzalez, A. Gibson, E. M. C. Hillman, R. M. Yusof, N. Everdell, D. T. Delpy, G. Zaccanti, and F. Martelli, “Assessment of an in situ temporal calibration method for time-resolved optical tomography,” J. Biomed. Opt. **8**(1), 87–92 (2003). [CrossRef]

**20. **R. M. Yusof, J. C. Hebden, A. Gibson, N. Everdell, T. Austin, J. H. Meek, S. R. Arridge, J. S. Wyatt, and D. T. Delpy, “Validation of the use of homogenous reference phantoms for optical tomography of the neonatal brain,” Proc. SPIE **4955**, 6–11 (2003).

**21. **T. Durduran, R. Choe, J. P. Culver, L. Zubkov, M. J. Holboke, J. Giammarco, B. Chance, and A. G. Yodh, “Bulk optical properties of healthy female breast tissue,” Phys. Med. Biol. **47**(16), 2847–2861 (2002). [CrossRef]

**22. **J. Selb, A. M. Dale, and D. A. Boas, “Linear 3D reconstruction of time-domain diffuse optical imaging differential data: improved depth localization and lateral resolution,” Opt. Express **15**(25), 16400–16412 (2007). [CrossRef]

**23. **M. Kacprzak, A. Liebert, P. Sawosz, N. Zolek, D. Milej, and R. Maniewski, “Time-resolved imaging of fluorescent inclusions in optically turbid medium – phantom study,” Opto-Electron. Rev. **18**(1), 37–47 (2010). [CrossRef]

**24. **A. Dalla Mora, A. Tosi, F. Zappa, S. Cova, D. Contini, A. Pifferi, L. Spinelli, A. Torricelli, and R. Cubeddu, “Fast-gated single-photon avalanche diode for wide dynamic range near infrared spectroscopy,” IEEE J. Sel. Top. Quantum Electron. **16**(4), 1023–1030 (2010). [CrossRef]