2.1.

Experimental Setup

Figure 1 reports the setup used for the experiments. The light source was a mode-locked laser (Fianium Ltd., United Kingdom) providing optical pulses [duration: 26-ps full-width at half maximum (FWHM)] at 820 nm at the repetition rate of 40 MHz. The light passed through a variable optical attenuator and was then injected into the sample through an optical fiber [core diameter: $200\mu \mathrm{m}$; numerical aperture (NA): 0.22; Lightech srl, Italy]. The photons re-emitted from the sample were then collected at two different points by means of two optical fibers (core diameter: 1 mm, NA: 0.37; Thorlabs Gmbh, Germany). We decided to use two detection points so as to increase the quality of the reconstructed image by acquiring data from multiple points of view. The distance between injection and collection fibers was set to 30 mm. The light exiting from each collection fiber was then imaged onto a SiPM detector by a set of two lenses (achromatic doublets lens, focal length: 30 mm) providing a $1\times$ magnification.

Fig. 1

Experimental setup: a pulsed laser light is attenuated and shined onto a phantom through an optical fiber. The backscattered photons are detected by the two SiPM modules, which provide the start signal to the TCSPC boards, while the stop one is given by the laser.

Each module embeds a $1{\mathrm{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. Additionally, the module hosts a thermoelectric cooler controller board, which is set to keep the detector at about 15°C. This value represents the best compromise to limit both primary dark count rate and afterpulsing (i.e., noise contribution due to charge carriers trapped during the avalanche ignition). For further details about the SiPM module design and characterization, see Ref. 27. Each SiPM module was then coupled to a TCSPC board (SPC-130, Becker&Hickl GmbH, Germany) to provide the start signal for time-stamp generation. The stop pulse was given to both TCSPC boards by the “synchronization” output of the laser.

To acquire the instrument response function (IRF) of the system, we directly faced the source fiber and one detector fiber at a time, putting a thin Teflon layer in-between to ensure the isotropic illumination of the latter. The IRFs obtained for a 10-s measurement time and after the subtraction of the mean value of the background noise for both SiPM modules are reported in Fig. 2. The peak timing jitter (FWHM) is around 260 ps for both detectors (262 ps for module 1 and 255 ps for module 2) while the tail with a long time constant ($\sim 3\mathrm{ns}$) starts 2.5 decades below the peak, in agreement with Ref. 32.

Fig. 2

IRF recorded with both SiPM modules (constant background subtracted). The region of interest (ROI) is the portion of the curve used for simulations.

The homogeneous liquid phantom was a water dilution of black Indian ink (which acts as an absorber) and Intralipid (providing scattering property). The proportion of the components were chosen so as to have an absorption coefficient (${\mu }_{\mathrm{a}}$) of $0.07{\mathrm{cm}}^{-1}$ and a reduced scattering one (${\mu }_{\mathrm{s}}^{\prime }$) of $12{\mathrm{cm}}^{-1}$ at 820 nm.

As an absorbing inhomogeneity within the homogeneous medium, we used five totally absorbing inclusions [black cylinders made of black polyvinyl chloride (PVC)] of different volumes. As reported in previous studies,^35,36 small totally absorbing objects properly reproduce realistic absorption inhomogeneity in a reliable and reproducible way and for a wide range of applicability. In particular, for a given ${\mu }_{\mathrm{s}}^{\prime }$ background, it is possible to identify a whole equivalence class of absorbing perturbations, with different size and shape, which yield the same temporal perturbation for any choice of geometry (e.g., reflectance, transmittance, source–detector distance), position of the object within the medium, and background absorption. Each equivalence class also contains a totally absorbing object whose volume [equivalent black volume (EBV)] can be taken as representative of the whole class. Some breaking of the equivalence are observed only when black object has a very small volume ($\le 50{\mathrm{mm}}^3$) and the object is quite close (depth $\le 10\mathrm{mm}$) either to the source or to the detector. The black inclusions we used in this work were cylinders with height equal to the diameter with values of 3.2, 4, 5, 6.8, 8.6 mm, a volume $\mathrm{EBV}=25$, 50, 100, 250, $500{\mathrm{mm}}^3$, and corresponding to an equivalent $\mathrm{\Delta }{\mu }_{\mathrm{a}}=0.056$, 0.087, 0.15, 0.37 and $0.94{\mathrm{cm}}^{-1}$ for a $1{\mathrm{cm}}^3$ volume, respectively.³⁶ In this paper, the inclusions will be named as “inclusion 1, 2, 3, 4, or 5” depending on their volume (e.g., we will refer to the $100{\mathrm{mm}}^3$ as “inclusion 3”). For what concerns data presented in the following section, in most cases, there is negligible discrepancy in the equivalence relation: only in two cases ($\mathrm{EBV}=25$, $50{\mathrm{mm}}^3$, with $\text{depth}=10\mathrm{mm}$), some 10% to 20% discrepancy in the amount of perturbation as compared to the equivalent realistic inhomogeneity is observed.³⁶

Since multiple measurements are needed to obtain a tomographic reconstruction during measurements, the injection and collection fibers [represented, respectively, as yellow star and open red circles in Fig. 3(a)] 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 plate in which holes for fiber tips were drilled. Using a motorized three-axes translation stage, the inclusion was moved over an area of $3.75\times 3{\mathrm{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\mathrm{mm}$ to 15 mm for $x$ and $y$ direction, respectively]. For each inclusion position [small black dots in Fig. 3(a)], we acquired the signal for 15 s (15 repetitions of 1 s each). As reference measurement, for each line, we moved the inclusion far enough ($x=-40\mathrm{mm}$) to consider its effect negligible.

Fig. 3

(a) Geometry of the scan in the $x-y$ plane: injection and collection fibers (yellow star and red open circle, respectively) are fixed while the inclusion (gray filled circle) was moved in the different positions (black dots). (b) Example of a 3-D ${\mu }_{\mathrm{a}}$reconstruction of the inclusion at 20-mm depth using simulation data.

For each perturbation, we first aligned the inclusion at 7.5 mm from the injection fiber in both $x$ and $y$ directions, and this point was set as origin for the scanning coordinates system. This choice was done to obtain a scanning area not centered to the inclusion position, thus being able to disentangle the reconstruction of the inclusion from possible artifacts arising in the center of the image.

Scans of all the inclusions in the $x-y$ plane were done at different depths (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. Also, only for “inclusion 3,” 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.

2.2.

Data Analysis

In the following, the diffusion approximation was used to model light propagation.³⁷ Data analysis was performed considering two sets of measurements $M_{sd}^A(t)$ and $M_{sd}^B(t)$, acquired on a reference medium $A$ (homogeneous medium of known absorption ${\mu }_{\mathrm{a}}^A$) and on the unknown medium $B$ (with unknown absorption ${\mu }_{\mathrm{a}}^B$, but estimated at iteration $k$ by ${\mu }_{\mathrm{a}}^{B(k)}$). As explained in Ref. 13, the reference measurement is needed to take into account the effect of the IRF. In the previous expression of measurements, $t$ is the time, $s$ is an index on the source position, and $d$ is an index on a detector position. As in Ref. 13, we used the combination $Y_{sd}(t)$ obtained with these measurements and known Green’s functions $G^A$ of medium $A$ and Green’s functions $G^{B(k)}$ of estimated medium $B$ at iteration $k$, which is linearly related to the absorption update $\delta {\mu }_{\mathrm{a}}^{(k)}={\tilde{\mu }}_{\mathrm{a}}^B-{\mu }_{\mathrm{a}}^{B(k)}$. Green’s functions $G$ are defined as solutions of the diffusion equation for a Dirac source. In the reported experiment, ${\mu }_{\mathrm{s}}^{\prime }$ was measured ($12{\mathrm{cm}}^{-1}$) and was not considered as an unknown parameter in the data analysis. After space discretization with finite volume method and time discretization with the MLTs,^13,38 we ended up with this system:

3.1.

Simulations

Figure 4 shows the reconstructed maps (in the $z-y$ and $x-y$ planes, first and second row, respectively) for the “inclusion 3” at depths ranging from 10 mm up to 30 mm (columns). The black circle reported in each map represents the volume of $1{\mathrm{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 $x-y$ slices at almost all depths, the reconstructed volume is larger than the expected $1{\mathrm{cm}}^3$ and, more precisely, it progressively increases at higher depths. This is due to the fact that the information of deep inclusions is carried by photons with a long time-of-flight, thus providing a lower spatial resolution.⁴¹

Fig. 4

Simulation: sections of 3-D reconstructions of a perturbation equivalent to “inclusion 3” at different depths.

With the aim to study the sensitivity to the amount of absorption perturbation, Fig. 5 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 1). In addition, we clearly see that as the real $\mathrm{\Delta }{\mu }_{\mathrm{a}}$ of the perturbation increases, the retrieved $\mathrm{\Delta }{\mu }_{\mathrm{a}}$ increases as well (even if the absolute value is lower than expected) while the reconstructed volume does not, as expected. Indeed, since the inclusions depth is the same, the spatial resolution of the system is expected to be the same.

Fig. 5

Simulation: sections of 3-D reconstructions for perturbations equivalent to the five inclusions at 15-mm depth.

From simulations, we can conclude that SiPMs are suitable for the tomographic reconstruction of even small absorption perturbations ($\mathrm{\Delta }{\mu }_{\mathrm{a}}=0.056{\mathrm{cm}}^{-1}$) buried down to 20 mm (see Table 1), with good accuracy in terms of depth and lateral resolution. Indeed, for both $x$ and $y$ localization, the recovered position always differs of $<0.5\mathrm{mm}$ from the expected one, whereas for the $z$ axis, the relative error in depth localization is always $<3\mathrm{mm}$ (for inclusion 2 at depth 15 mm, which corresponds to an error of 16%). By contrast, the absolute quantification of the perturbation is not accurate, and its accuracy decreases upon the increase of the inclusion’s volume. However, the simulations done with the quasi-ideal detector reported in Table 1 clearly demonstrate that this lack in accuracy is not dependent on the peculiar features of the SiPM (e.g., slow diffusion tail), since similar results in quantification are recovered also using a state-of-the-art detector.

Due to the encouraging results obtained in the simulations, we have studied the suitability of the SiPM for TD-DOT with experimental measurements.

3.2.

Experimental Measurements

The effect of the depth (ranging from 10 to 30 mm, columns) for “inclusion 3” in phantom measurements can be seen in Fig. 6, where the $z-y$ (first row) and $x-y$ (second row) maps are reported. Up to 25 mm depth, the spatial localization (depth and $x-y$) is comparable with simulations. 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 than in simulations, but this is most likely due to instabilities of the experimental system (e.g., laser instability, ageing 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{\mathrm{cm}}^{-1}$ (for $1{\mathrm{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.

Fig. 6

Experiment: sections of 3-D reconstructions of the “inclusion 3” at different depths.

To ascertain the sensitivity on the amount of absorption perturbation, Fig. 7 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 terms of depth and $x-y$ position.

Fig. 7

Experiment: sections of 3-D reconstructions for the five inclusions at 15-mm depth.

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 $\mathrm{\Delta }{\mu }_{\mathrm{a}}$ upon an increase of the depth.

As can be seen in Table 1, the depth localization is quite accurate for all depths with an error always lower than 2.5 mm (up to 20-mm depth) that corresponds to a maximum error of 17%. We can also notice that experiments and simulations show a good agreement with a discrepancy in the position of the inclusion, which is always lower than 2 mm for depths (except for 30 mm).

The retrieved coordinates of the $x$ position of the inclusion recovered from the experiments seem 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\mathrm{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 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 on the retrieved $\mathrm{\Delta }{\mu }_{\mathrm{a}}$, measurements and simulations are usually in good agreement (see Table 1). The only exception is for inclusion 1 and 2 at a shallower depth, which could be partially ascribed to the deviation of the equivalence relations between EBV and the realistic $\mathrm{\Delta }{\mu }_{\mathrm{a}}$, which, for these particular cases, implies $\sim 20\%$ higher contrast for the black object in the experiments. On the other hand, there is a clear underestimation in the reconstructed $\mathrm{\Delta }{\mu }_{\mathrm{a}}$ upon increasing the depth. Since this trend 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. Indeed, as discussed in Sec. 3.1, simulations using an almost ideal detector (without long tail) are substantially in agreement with SiPM simulations. Thus, we can conclude that SiPM performances do not alter $\mathrm{\Delta }{\mu }_{\mathrm{a}}$ quantification capabilities as compared to a state-of-the-art detector.

The absolute quantification of the perturbation is a well-known problem in tomography⁴² in particular for deep inhomogeneities. More comprehensive studies are needed to quantify this effect properly, identify the determinants, and try to improve quantitation. Surely, the decrease in retrieved ${\mu }_a$ for increasing depth is a concern when imaging complex structures since fluctuations in reconstruction due to shallower inhomogeneities can hinder deeper and attenuated structures. From the data in Table 1, this effect—for the background properties presented in this study (${\mu }_a=0.07{\mathrm{cm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=12{\mathrm{cm}}^{-1}$)—is around a factor of 3 to 4 between a depth of 10 and 20 mm. Yet, this effect is compensated by the broadening of the reconstructed inhomogeneity, which implies that the overall impact of the deeper perturbation on the image is somehow augmented. In a different study, with conventional detectors, we have shown that fixing the size of the reconstructed inhomogeneity to the true value, as could be obtained through a coregistration with another imaging modality, completely solves the depth reduction in reconstructed ${\mu }_{\mathrm{a}}$.⁴³