2.1.

Setup

The imaging setup consisted of a $512\times 512\text{pixel}$ SPAD array [Fig. 2(a)] (SwissSPAD2, EPFL, Neuchâtel, Switzerland) with 10.5% fill factor, coupled to a 635 nm 20 MHz pulsed laser (LDH-D-C-635M laser diode and PDL 828 driver, Picoquant, Germany) delivering ultrashort laser pulses ($<120\mathrm{ps}$). The main sync signal generated by the laser driver was converted to a +5 V TTL signal provided as an input to the camera’s field-programmable gate array (FPGA) (Opal Kelly XEM7360, Xilinx, San Jose, California, United States). The laser was operated at maximum intensity, providing an optical power density of $160\mu \mathrm{W}/{\mathrm{cm}}^2$ in the target plane. The SPAD device and its FPGA were mounted in a custom aluminum case with forced-air cooling. A C-Mount was attached to the case and aligned with the SPAD sensor to relieve stress on the sensor’s printed circuit board. A 50 mm, 1:1.8 D lens (Nikkor, Nikon, Tokyo, Japan) was used in a first experiment and a 50 mm F0.95 lens (TV-lens, Lenzar, Jupiter, Florida, United States) in a second experiment to increase light collection capabilities of the system.

Fig. 2

(a) Layout of the SwissSPAD2 sensor. (b) Time diagram of a laser pulse’s TP acquisition sequence. Both the laser and the camera are triggered by a 20 MHz clock signal. The number of gate sequences is defined by the user, and each of these sequences can be composed of a chosen number of frames (not shown here for clarity). (c) A representative curve depicting the SPAD sensor’s linear range for fixed acquisition parameters.

For each experiment, the SPAD camera was placed directly above the target plane at a working distance of 50 cm. The laser beam was provided as input to the imaging setup using a fiber with a lens and a top hat diffuser fixed on the fiber terminus. The laser was placed in a way to reduce its angle with the normal to the sample’s surface while ensuring that it did not obstruct the camera’s field of view. Data transfer was achieved via USB from the FPGA to a dedicated imaging system computer.

2.1.2.

Data acquisition

The time-resolved acquisition sequence described above results in a 3D image stack with dimensions of $256\times 472\times 1400$, corresponding to the spatial rows, spatial columns, and time points, respectively. These data were then analyzed by selecting spatially averaged regions of interest (ROIs) with constant sizes, across the entire time scan, resulting in $1\times 1400$ TPs. The following section details how these TPs were acquired and processed.

Time profile generation

Depending on the signal-to-noise ratio (SNR) and lateral spatial resolution desired, TPs were generated over varying ROI sizes. Increasing the TP ROI size increases SNR at the cost of lower spatial resolution and vice versa. For each pixel, a TP corresponded to the number of photons detected as a function of time. Each TP was normalized to its maximum intensity, corresponding to the average value of the 200 most intense frames, and to its minimum intensity, corresponding to the average value of the 200 last frames where no fluorescent signal was detected; this normalization mapped TP values between the average maximum and minimum to the range of [0, 1] [see Fig. 3(a)].

Fig. 3

(a) Fluorescence TPs from inclusions in 3% Intralipid solution at depths of 0, 2, and 4 mm, normalized to maximum intensity for comparison. For each plot, the SNR is displayed for an ROI corresponding to an $8\times 8$ binned pixel area. (b) Corresponding intensity images showing the spatial blur and signal loss with increasing inclusion depth. (c) From top to bottom, IRF from white paper, fluorescence intensity image, and reflectance intensity image averaged over 100 frames during the pulse; the co-registered ROIs enable fluorescence correction. (d) How fluorescence TPs are corrected by corresponding reflectance TPs. Setup changes induce time shifts in the TPs, and correcting for this allows accurate comparison of TPs between acquisitions. IRF, instrument response function.

Instrument response function

The IRF of the SPAD camera provided a reference for data fitting and correction in each experiment. The IRF served to rectify pixel temporal inconsistencies as well as variations in illumination across the field of view. Each fluorescence or reflectance data TP fitting was derived using a corresponding ROI on the IRF time scan [see Fig. 3(c)]. In each setup configuration, the IRF was acquired to account for the setup geometry, because any adjustment in setup geometry modified the IRF both spatially, because of the laser angle, and in time, because of the camera-to-sample working distance. The IRF was acquired using a sheet of white paper placed in the imaging plane of the setup. Here this paper sheet was considered a Lambertian surface, allowing for temporal correction of the fluorescence and reflectance signals over the entire field of view.

Fluorescence data

IRDye 680RD was excited using a 635 nm wavelength laser, corresponding to the dye’s first absorption Q-band. A 697 nm band-pass filter with 75 nm bandwidth and $\mathrm{OD}>6$ (Edmund Optics, Barrington, New Jersey, United States) separated the emission signal from the excitation wavelength.

Reflectance data

Reflectance data were acquired in the same way as fluorescence data, except that the 697 nm band-pass filter was removed, and linear polarizing filters were placed on the laser fiber optic output and on the front of the camera lens. The polarizers were offset by 90 deg. It is well-known that polarized light remains polarized after undergoing specular reflection, whereas polarized light loses its polarization when passing through turbid media. Therefore, the use of crossed polarizers minimized specular reflections from the surface of each sample. This was especially important when working with surfaces that exhibited high levels of specular reflection, such as solutions and gelatin phantoms, as these reflections can otherwise overpower the useful signal originating from backscattered light. Additionally, the polarizers functioned as ND filters to account for the absence of the spectral filter and to prevent sensor saturation.

Acquisition sequence

For each imaging session, an IRF was acquired. All data collected during the session were corrected using the same IRF. To compensate for samples’ surface-to-sensor distance changes and surface geometry, each fluorescence TP acquisition was coupled with a reflectance TP acquisition using the same settings. To facilitate data collection for untrained operators, a graphical user interface (GUI) was implemented using MATLAB (vR2018a, MathWorks Inc., Natick, Massachusetts, United States). The GUI allowed for automatic acquisition of a set of five fluorescence TPs and five reflectance TPs, with a break in between allowing the user to switch the filters. Raw data were automatically saved for subsequent analysis. An entire acquisition sequence using the GUI took $<10\mathrm{s}$.

2.2.

Sample Preparation

2.3.

Data Processing

2.3.1.

Time profile fitting

To obtain accurate TOF measurements, it was necessary to fit both fluorescence and reflectance TPs to an analytical model. The analytical model was derived by examining the various stages of how light pulses interact with the medium and different layers of the geometry. The IRF was characterized as the TP resulting from a laser pulse propagating through a negligible thickness of nonscattering medium (i.e., air), reflecting off a Lambertian surface, and subsequently propagating back to the detector. In the case of a laser pulse propagating through a turbid medium, the TP is altered, and a term representing the medium transfer function must be incorporated. This is expressed by the following equation:

Eq. (1)

$$\mathrm{TP}=\mathrm{IRF}*H_{\text{medium}}{*\mathrm{\Delta }}_{\text{air}},$$

where TP is the resulting time profile registered by the sensor, IRF is the instrument response function described above, $H_{\text{medium}}$ is the medium transfer function, ${\mathrm{\Delta }}_{\text{air}}$ is the delay corresponding to light’s travel through air [see Fig. 4(a)], from the laser head to the tissues and back (negligible), and * is the convolution operator.

Fig. 4

(a) Schematics of the setup with wavefronts represented by dotted lines to delimitate TOF paths. (b) Flowchart of depth calibration using TOF measurements and either reflectance measurements or SFDI measurements to correct for tissue optical properties.

If the pulse train, in addition to traveling through a turbid medium, encounters a Stokes shift due to a fluorescent emitter, another convolution must account for the phenomenon:

Eq. (2)

$$\mathrm{TP}={\mathrm{IRF}*\mathrm{H}}_{\text{medium}}*e^{-t/{\tau }_{\mathrm{fluo}}}{*\mathrm{\Delta }}_{\text{air}},$$

where $\tau$ is the fluorescent lifetime of the probe assuming a monoexponential response.

Reflectance data

Reflectance data can be fitted using Eq. (1), because reflectance results primarily from light interactions with the superficial layer of a turbid medium,³² and it can be assumed that the fluorescent contribution is negligible in terms of intensity. In this study, we endeavored to identify a mathematically simple yet effective expression that could enhance computational efficiency. The medium transfer function $H_{\text{medium}}$ can be effectively modeled through a Gaussian function, as empirically substantiated by the reflectance fitting results. The mean goodness of fit for reflectance data was found to be $0.997\pm 0.011$, providing strong evidence for the appropriateness of the Gaussian function in capturing the characteristics of $H_{\text{medium}}$. The TP fitting model can therefore be described as follows:

Eq. (3)

$$I_{\mathrm{fit}}(t)=\mathrm{IRF}*\frac{A}{\sigma \sqrt{2\pi }}{e}^{\frac{-{(t-{\mathrm{\Delta }}_{\text{tissue}})}^2}{2{\sigma }^2}},$$

where $I_{\text{fit}}$ is the analytical model fitting the reflectance TP, dependent on time $t$, and $A$, ${\mathrm{\Delta }}_{\text{tissue}}$, and $\sigma$ are the amplitude, mean, and standard deviation, respectively, of the Gaussian distribution modeling the medium transfer function. In practice, ${\mathrm{\Delta }}_{\text{air}}$ is corrected by the IRF. The three parameters $A$, ${\mathrm{\Delta }}_{\text{tissue}}$, and $\sigma$ were set as variables, and $I_{\text{fit}}$ was fitted to the raw data using a least squares method.

Fluorescence data

Following similar logic, the backscattered fluorescence TP fitting can be modeled using Eqs. (2) and (3), where the fluorescent term can be described as follows:

Eq. (4)

$$I_{\mathrm{fit}}(t)=\mathrm{IRF}*\frac{A}{\sigma \sqrt{2\pi }}{e}^{\frac{-(t-{\mathrm{\Delta }}_{\mathrm{tissue}}{)}^2}{2{\sigma }^2}}*e^{-t/\tau },$$

where $\tau$ is the fluorescence lifetime of the probe. For all the fitting done in this paper, $\tau$ was set to 0.7 ns, as measured using time correlated single-photon counting lifetime measurements (Fluoromax, Kyoto, Japan). This simple model is illustrated in Fig. 5(d).

Fig. 5

(a) Influence of depth on the fluorescence TP’s shift and dispersion. The figure shows the fluorescent profiles at three different depths, as well as the reflectance TP for a $1\times 1\times 4{\mathrm{cm}}^3$ IRDye 680RD carboxylate inclusion embedded in a 3% IL aqueous solution. (b) Influence of the medium’s scattering coefficient on the fluorescence TP shift and dispersion. The influence on the reflectance TP is less pronounced but can be seen in the inset window. (c) Schematic of the setup showing the fluorescent cuvette immersed in an Intralipid solution. (d) The TPs analytical fitting model is the result of the convolution between the IRF, an exponential function corresponding to fluorescence decay, and a Gaussian function that reflects the medium’s transfer function. (e) Relative TOF ${\mathrm{\Delta }}_{\text{tissue}}$ compared to inclusion depth at varying IL concentration. ${\mathrm{\Delta }}_{\text{tissue}}$ corresponds to the Gaussian’s central value displayed in (d). The red plot represents the linear relation between the slope of each of the five TOF curves as a function of IL concentration. IL, Intralipid.

2.3.2.

Time profiles comparison

To compare results from different acquisitions—i.e., experiments where modifications were made to the setup or where samples were replaced (as was the case for the liquid phantom and the gelatin slabs experiments), it was necessary to ensure that any changes did not introduce bias. Thus to accurately compare fluorescent TPs from different acquisitions, each fluorescent TP was time-corrected to its corresponding reflectance TP’s rising edge. Raw reflectance TP rising edges were fitted with linear functions, as shown in Fig. 3(d), which were then used as a time reference. Figure 3(d) provides an illustration of the fluorescence TP before and after correction. The red lines represent the linear fit of the reflectance rising edge, and $\mathrm{\Delta }t$ represents the time delay applied to correct the fluorescence TPs. The corresponding IRF is shown for reference. To enhance clarity, all TP data shown in Figs. 3, 5, and 6 were smoothed with a 35 element (i.e., 623 ps) moving average filter.

Fig. 6

(a) Reflectance TPs with respect to IL concentration. Although temporal variations are small relative to the overall pulse length, the inset window highlights how IL concentration impacts the falling edge of each pulse. (b) SFDI-derived reduced scattering coefficient spectra for different IL concentrations. (c) Illustration of the linear relation between measurements from (a) and (b) and IL concentration. ${\mathrm{\Delta }}_{\mathrm{refl}}$ represents the time from reflectance TPs at 2% of the normalized pulse. (d) Agreement is observed between SFDI-based and reflectance-based depth measurement calibration.

2.4.

Depth Calibration/Correction

It is important to note that an absolute TOF measurement may not necessarily reflect an accurate distance traveled by photons. This is because the TOF of backscattered photons is influenced by the optical properties of the turbid medium. To obtain precise depth measurements, it is essential to calibrate the absolute TOF measurements in accordance with the medium’s optical scattering and absorption properties. The two methods described below were evaluated for the purpose of performing this calibration and are described in Fig. 4(b).

3.1.

Liquid Phantom

The depth sensing system was first tested on a simplified tissue phantom model. The phantom consisted of a cuvette immersed in turbid solutions with varied scattering coefficients as shown in Fig. 5(c). TPs shown in Figs. 5(a) and 5(b) were obtained as described in Sec. 2.3.1 using a $40\times 40\text{pixel}$ ROI. Figure 5(a) shows the influence of the inclusion depth on the fluorescence TPs recorded by the SPAD camera. These pulses are plotted next to corresponding IRF and reflectance TPs. Because the reflectance TPs were very similar for the various inclusion depths, only one reflectance curve is displayed in Fig. 5(a) for clarity. The plot shows that the fluorescence TP is affected in two ways as the inclusion depth increases: first, the TPs shift toward the right as depth increases, i.e., ${\mathrm{\Delta }}_{\text{tissue}}$ increases, and second, the TP becomes increasingly dispersed with depth, i.e., $\sigma$ increases. The pulse spread is not clearly visible in these figures but will be described in more detail in Sec. 3.5. Figure 5(b) shows the influence of the medium’s optical scattering level at 5 mm depth, expressed here in terms of IL concentration, IL%, on the pulse. The impact of IL% on the pulse is similar to that observed in Fig. 5(a): an increase in scattering coefficient leads to greater temporal shifts and dispersion of the pulse. The coefficient ${\mathrm{\Delta }}_{\text{tissue}}$ expressed in Figs. 5(a) and 4(b) corresponds to the Gaussian’s mean value fitting coefficient. To aid the comprehension of fluorescent pulse fitting, an illustration of the convolution product is presented in Fig. 5(d). The model used to fit raw data results from the convolution of the IRF, an exponential decay representing the fluorescent re-emission, and a Gaussian function describing the medium’s effect on the light pulse TP. Later in this paper, we will demonstrate that ${\mathrm{\Delta }}_{\text{tissue}}$ is the most accurate reporter of TOF (see Sec. 3.5). Figure 5(e) finally summarizes the results obtained from five different datasets, each corresponding to an IL% ranging from 0.5% to 3.0%. TOF is expressed as ${\mathrm{\Delta }}_{\text{tissue}}$ and is compared to the inclusion depth. A linear relationship can be observed between TOF and depth for a given IL%. Once fitted, the slopes $k_{\mathrm{fluo}}$ of these five curves can be plotted against IL%, revealing a linear trend.

3.2.

System Calibration

This section focuses on absolute TOF measurement calibration by considering the optical properties of the bulk medium. Figure 6 illustrates the results of TOF calibration using two different methods for assessing the optical properties, as described above. The first method [Fig. 6(a)] utilizes reflectance TP measurements. Although reflectance TPs are less affected by tissue optical properties than fluorescence TPs, some variation can still be observed for varying IL concentrations. Figure 6(a) shows reflectance TPs obtained from IL in aqueous solution with 1%–5% IL concentrations, along with the corresponding IRF. All five reflectance TPs were aligned on their rising edge, as described in Sec. 2.3.2, and compared at a 2% threshold of the normalized pulse, as described in Sec. 2.4.1. The corresponding time value is denoted as ${\mathrm{\Delta }}_{\mathrm{refl}}$, which exhibits a linear relationship when plotted as a function of IL%. A closer look is provided in the inset window in Fig. 6(a). The portion of this curve used for calibration is shown in Fig. 6(c) as the blue line. It should be noted that ${\mathrm{\Delta }}_{\mathrm{refl}}$ is a relative time value, where ${\mathrm{\Delta }}_{\mathrm{refl}}=0$ was chosen arbitrarily at $\mathrm{IL}\%=0.5$.

The second method used for characterizing the medium’s optical properties was SFDI. Figure 6(b) shows the reduced scattering coefficient ${\mu }_s^{\prime }$ corresponding to various IL solutions in aqueous solution, measured with the SFDI imaging system. ${\mu }_s^{\prime }$ at 659 nm is plotted as a function of IL% in Fig. 6(c), showing only the range of IL% used in the liquid phantom experiment (Sec. 2). The linear relationship between ${\mu }_s^{\prime }$ and IL% is evident.

Finally, Fig. 6(d) shows the data obtained from the IL solution phantom presented in Fig. 6(e), corrected using the reflectance data (blue curve) and SFDI-based measurements (green curve), respectively. Both lines represent the trendline of the average corrected data, and the error bars represent the standard deviation of the five measurements. These data are corrected using Eqs. (10) and (11), respectively.

3.3.

Gelatin Slabs Phantom

To create a more controlled and realistic scenario for imaging fluorescent inclusions at varied depths, we conducted a third experiment using gelatin slabs that were 1 mm thick and had tissue-like optical properties. These slabs were stacked on top of each other, as illustrated in Fig. 7(a). The phantom setup allowed for precise control of the inclusion depth by adding or removing layers one by one. The schematic from Fig. 7(a) shows two scenarios, one representing a 2 mm thick layer above the inclusion (${\mathrm{\Delta }}_{\text{tissue}1}$) and the second representing a 5 mm thick layer above the inclusion (${\mathrm{\Delta }}_{\text{tissue}2}$). The arrows illustrate the path of a photon through the tissue, toward the dyed phantom portion, and back to the surface. In the same way as Fig. 6(e), Fig. 7(b) shows TOF as a function of depth, for different IL concentrations. In this case too, a linear relationship was found between TOF and inclusion depth for each IL concentration. These data were corrected for the medium’s optical properties using only SFDI-derived optical property measurements.

Fig. 7

(a) Illustration of the gelatin slabs phantom, including a drawing that shows the layering of the dyed and undyed slabs as well as the time difference between photons traveling through different thicknesses of the medium. The photograph illustrates the process of peeling away slabs to adjust the inclusion depth. (b) Relative TOF ${\mathrm{\Delta }}_{\mathrm{tissue}}$ versus inclusion depth at varying IL concentrations from the gelatin slabs phantom. (c) Depth correction obtained using SFDI data from the varying IL phantoms.

3.4.

Depth Correction

Figure 8 shows the result of surface luminance correction for depth attenuation in both the liquid phantom experiment [Fig. 8(a)] and the porcine gelatin slab phantom [Fig. 8(b)]. In both graphs, fluorescent intensity was measured in terms of photon counts at the phantom surface for varying inclusion depths. These intensities were then normalized to the max intensity for each IL concentration and plotted as solid lines. These data were then corrected using Eq. (10) and plotted as the dotted lines in both graphs. Figure 8(b) demonstrates more consistent intensities as the inclusion depth increases, as a result of depth attenuation correction.

Fig. 8

Surface fluorescence intensity at varying inclusion depths (full lines) and corrected intensity representing the luminance at depth (dotted lines) for varying IL concentrations for inclusions embedded in (a) IL aqueous solution and (b) gelatin/blood/IL mixture.

3.5.

Plastic Tube in Gelatin Phantom

The last experiment was intended to realistically simulate a fluorescent inclusion in tissue to test the depth calibration developed in the previous section. For this purpose, two fluorescent cylinders were randomly inserted at different depths into a tissue-like phantom composed of porcine gelatin, blood, and IL (see Sec. 2.2.2). The results from the TOF imaging of this phantom are shown in Fig. 9. Figure 9(d) shows a cross section of the phantom. On the bottom is a white light image of the phantom cut in half after imaging, and on the top is a fluorescence image of the same cross section, captured with the Odyssey CLx imaging system. (The white light image does not correspond to the cross-section imaged in this experiment and is only shown for illustration purposes.) Although air gaps can be seen between the inclusion and the surrounding bulk material in the cross-sectional image [Fig. 9(d)], these gaps are a result of cutting the phantom for the representative white light image. Importantly, these gaps were not present during the TOF imaging. Figures 9(a)–9(c) show the results of the measured TPs fitting parameter $\sigma$ [Fig. 9(a)] and ${\mathrm{\Delta }}_{\text{tissue}}$ [Fig. 9(b)] using $8\times 8\text{pixel}$ binning. In this case, the $8\times 8\text{pixel}$ binning reduced noise and increased contrast. The SNRs for pixels corresponding to the shallowest and deepest inclusions were $\sim 10$ and $\sim 2$, respectively, in the case of single pixels, and $\sim 74$ and $\sim 19$ when employing the $8\times 8\text{pixel}$ binning strategy. The Gaussian fitting parameter $\sigma$ is a good reporter of the inclusion $x-y$ position in the target plane [Fig. 9(a)]. The inclusion closer to the surface shows up distinctly while the deeper inclusion is visible but less pronounced. Figure 9(b) shows the TOF map corresponding to the fitting factor ${\mathrm{\Delta }}_{\text{tissue}}$. It appears that ${\mathrm{\Delta }}_{\text{tissue}}$ is a good reporter of TOF and shows inclusion depth in more detail. The time scale used for the color bar represents a time relative to the first returning light pulse. Thus the area corresponding to the shortest distance between the inclusion and the phantom surface corresponds to ${\mathrm{\Delta }}_{\text{tisue}}=0$. On the left and right side of the TOF images, the areas showing high noise correspond to the limit of the phantom bulk. Since these areas do not exhibit fluorescence, the fitting model provides inconsistent results as shown on the pixel goodness-of-fit map [Fig. 9(c)]. To overcome this problem, the $8\times 8$ binned TOF image was smoothed using a spatial Gaussian filter [Fig. 9(e)] and overlayed on top of the fluorescence intensity image [Fig. 9(f)]. The goodness-of-fit consistently exhibited values of $0.99\pm 0.03$ for pixels corresponding to phantom areas, with these values representing the mean and standard deviation of the $R^2$ values, respectively. The resulting image [Fig. 9(f)] shows the final TOF image free of background noise. A 3D representation of the normalized depth is shown in Fig. 9(g), where TOF has been translated to a relative depth from the phantom surface. Figure 9(h) shows the topography obtained from the TOF measurement and corrected using the slab phantom calibration described in Sec. 3.3.

Fig. 9

Fluorescent inclusions embedded in the tissue-like gelatine phantom with cylindrical inclusions. (a) Map of Gaussian fitting parameter $\sigma$ giving a qualitative idea of how much signal coming from two fluorescent tubes gets blurred by the medium. (b) Map of fitting coefficient ${\mathrm{\Delta }}_{\mathrm{tissue}}$, which is a surrogate measurement of TOF. (c) Goodness-of-fit for each $8\times 8$ binned pixel expressed as $R^2$. (d) White light image and corresponding fluorescence image of the gelatine phantom cross section. (e) Averaged fluorescence intensity image with histogram equalization for visibility. (f) Overlay of the $8\times 8$ binned TOF map on the equalized intensity picture, which provides background noise cancellation. (g) Topological representation of overlay from (f). (h) Final measured topology compared to true topology overlayed on the phantom cross section.

3.6.

Tumor Resections

Figure 10 shows the result of two labeled specimens from two freshly resected head and neck tumors (details in Sec. 2.2.3). Figure 10(a) shows fluorescence intensity overlayed on a reflectance image of the two tumor slices. The hot spots suggest inhomogeneous distribution of the imaging agent. Figure 10(b) shows the fluorescence TOF map expressed as the fitting factor ${\mathrm{\Delta }}_{\text{tissue}}$. The bright spots observed in the background reflectance image are due to specular reflections. (Cross-polarization filtering was not implemented for this experiment.) Figure 10(c) shows the fitting factor $\sigma$ representing the excitation pulse’s time spread.

Fig. 10

(a) Fluorescent signal from two representative breadloaf slices from two freshly resected head and neck tumors expressed as SNR from an $8\times 8$ binned pixel image, overlayed on top of the reflectance image. (b) Relative fluorescence TOF expressed as ${\mathrm{\Delta }}_{\text{tissue}}$. (c) Map of the fitting parameter $\sigma$ representing the light pulse time-spread.

4.1.

Data Processing

Fluorescence TOF depth sensing using a large-format SPAD camera involved multiple technical aspects, each of which is discussed in the following sections.

4.2.

Liquid Phantoms

The first part of this study aimed to determine if TOF as measured by the system was a good reporter of fluorescent inclusion depth in turbid media. The first set of experiments focused on measuring fluorescent TP from a phantom where the inclusion depth was controlled precisely. The liquid phantom consisted of IL in aqueous solution at varied concentrations. This phantom was a simplified case since its absorption was close to 0. That allowed for increased signal and improved TOF sensing at depth.

Even though inclusion depth was controlled precisely, Fig. 5(e) shows that the five linear plots do not converge at a depth of 0, which corresponds to the case where the inclusion was not surrounded by any liquid. The main reason for this discrepancy is that the first addition of surrounding liquid was not precise due to superficial tension at the surface of the sample. As a result, the depths shown in Fig. 5(e) are not completely accurate. However, once the superficial tension was broken, a controlled thickness of solution could be added, making the change in depth precise. For this reason, in Fig. 5(e), linearity is conserved from 1 to 5 mm, but a slight offset is observed between 0 and 1 mm. As a result, the data for depths between 0 and 1 mm were not included in the analysis. The coefficient $c_{\mathrm{fluo}}$ from Eqs. (10) and (11) is in theory constant and only depends on the setup geometry.

Another element that was not considered in this set of experiments is the 1 mm thickness of the plastic cuvette’s wall. Due to the speed of light in a transparent medium, it is not possible, given the temporal resolution available, to resolve such thin layers. Therefore, it was assumed that the plastic walls had a negligible influence on the time-resolved data.

The goal of this first experiment was to highlight the relationship between measured TOF and inclusion depth. It was shown that this relationship is linear in the observed depth range (1 to 5 mm). It also showed that linearity was maintained at various levels of optical scattering. The higher the scattering coefficient is, the faster TOF increases with depth is. This is shown by the red curve in Fig. 5(e), which reports the slope $k_{\mathrm{fluo}}$ as a function of IL concentration. One final take away from this first set of experiments was that the system’s ability to resolve sub-millimeter depth was due to the optical properties of the target medium (i.e., high scattering). The depth sensitivity achieved in this set of experiments would not have been possible otherwise, even considering the high refractive index of the medium.

4.3.

Depth Calibration

Depth calibration can be achieved by quantifying changes induced by varied optical properties, as described by Eqs. (10) and (11). Two methods were investigated for this purpose: first, reflectance-based calibration using direct time-resolved measurements from the SPAD camera; and second, using an established external modality (i.e., SFDI). Both calibration models were developed empirically by observing the relationship between different parameters. It should be noted that while these methods provide quantification of bulk optical properties in surface tissues, they do not possess the capacity to evaluate these parameters in terms of depth. Consequently, the general approach presented in this study operates under the assumption that variations in optical properties are insignificant along the $z$ axis.

4.4.

Gelatin Slabs Phantom

The second experiment involved a phantom with more realistic parameters compared to the first liquid phantom experiment. Porcine gelatin, bovine blood, and IL were used to fabricate a tissue-like phantom. The fluorescent inclusions were dyed with IRDye 680LT at a concentration of 50 nM (versus $10\mu \mathrm{M}$ for the liquid phantom). The 1 mm gelatin slabs offered precise control of the inclusion depth from 0 to 5 mm depth. Since the dye concentration was relatively low, the system was upgraded with a larger aperture lens and a shorter working distance. Results from Fig. 7(b) show that in the same way as in Fig. 5(e), there is a linear relationship between TOF and fluorescent inclusion depth. Data from this experiment were noisier than the data from the previous experiment. TOF data at 5 mm were not shown here, because signal reached the noise level near a depth of 4 mm due to the medium’s relatively high absorbance.

This experiment showed the limit of the SPAD LiDAR system. Even though TOF reconstruction is theoretically insensitive to fluorescence intensity, a low SNR is correlated with TP fitting inaccuracies. Averaging signal over a larger pixel area allowed for better results but at the cost of reduced spatial resolution.

4.5.

Plastic Tube in GELATIN Phantom

The third experiment involved a fluorescent tube embedded in gelatin. The phantom was fabricated in a way that the inclusion position in space was not precisely known. Therefore, this experiment more closely emulated a real-life scenario where the inverse problem had to be solved using the calibration based on the gelatin slab phantom. Upon analyzing the results depicted in Figs. 5(a) and 5(b), it appeared that the spread of TPs did not change with inclusion depth or IL concentration. However, a closer examination in Fig. 9(a) revealed that the pulse width, denoted as $\sigma$, was indeed affected by depth. The result was noisy, and the deeper inclusion was only slightly discernable in Fig. 9(a). TOF, expressed as ${\mathrm{\Delta }}_{\text{tissue}}$, on the other hand, seemed to be a more effective reporter of the inclusion position. When comparing Figs. 9(g) and 9(h), the relative distance between objects was conserved. For both fitting parameters, $8\times 8$ binned pixel images were used, as it improved contrast, particularly in the case of $\sigma$. For ${\mathrm{\Delta }}_{\text{tissue}}$, pixel binning did not have a significant impact on the result. The outgrowth seen at the bottom right of the shallowest tube was due to the transparent glue used to seal the tube. As it is transparent, the measured TOF is reduced as it acts like a waveguide; it does not represent a fluorescent inclusion.

To differentiate background noise from actual signal, the TOF map was Gaussian filtered and overlayed on the fluorescence intensity image. To ensure accuracy, only TOF measurements from regions expressing fluorescence were retained. However, this method has the disadvantage of discarding information in regions with low fluorescence intensity, such as the bottom right side of Fig. 9(f). This experiment aimed to simulate real condition measures where the spatial positions of inclusions were unknown, the goal being to locate them in space.

These combined results helped to verify that the methodology employed in this study would be robust for recovering the spatial positions of inclusions within tissue. However, the depth measurement capability has limits in cases where inclusions are located close to each other, due to scattering. Further deconvolution to account for high scattering is needed to obtain more accurate depth measurements.

4.6.

Tumor Resections

Finally, we were able to acquire SPAD LiDAR images of representative tumor slices from two freshly resected head and neck tumors, carried out as part of a larger ongoing study. In this case, signal levels were low, because the tumor was labeled with ABY-029 at lower (near microdose-level) concentration (SNR of $\sim 1.4$ for single pixels and SNR of $\sim 3$ for $8\times 8$ binned pixels). Furthermore, ABY-029 is an Affibody molecule labeled with IRDye 800CW, whose excitation and collection spectra did not align optimally with the parameters used in our system. Nevertheless, it was an excellent opportunity to image human tissue that had just been resected, albeit without being able to ultimately compare to the ground truth of depth from which the fluorescence came.

The $\sigma$ map [Fig. 10(c)] did not appear to be correlated with the TOF map. This could be due to both the low SNR of the dataset, providing poor fitting of fine parameters, and/or inhomogeneities in the tissue optical properties.

The tumor tissue samples provided a first assessment of how the SPAD LiDAR TOF approach performs when imaging biological samples with irregular surface profiles. While the true depths remain unknown, the system demonstrated consistent results in terms of TOF range and tumor distribution across the resection samples.

To address the challenges arising from the dataset’s low fluorescence intensities and pronounced noise levels, we implemented an 8 × 8 binning technique for the sensor. Although our initial aim in employing a large-format SPAD was to enhance spatial resolution, our findings revealed that pixel binning proved effective in overcoming these limitations and restoring valuable SNR. Through the consolidation of pixels via binning, we still achieved a sufficient spatial resolution, enabling the visualization of intricate biological details with a fine resolution down to 1 mm in scale—particularly evident in the context of tumor resection imaging.

4.7.

System Improvement and Opening

The use of a high-resolution, time-resolved SPAD camera as a LiDAR system was overall a success. Absolute time-resolved reading from the system was not quantitatively informative of true depth. To achieve true depth measurements, the system needs to be calibrated for the target medium’s specific distributions of optical properties. To be able to efficiently work this way, a lookup table with calibration corresponding to various combinations of optical properties would be needed. Machine learning could help address this aspect of system improvement.^38,39

In any case, the SPAD LiDAR method still requires a way to quantify medium’s optical properties. It was found that using SFDI was a robust way to do so, at the expense of requiring a separate imaging setup and additional processing time. Alternatively, time-resolved, patterned light projection composed of an array of dots could be used to recover medium’s properties.⁴⁰ This would offer a compact method to correct TOF data that leverages the same SPAD sensor.

Finally, increasing excitation power could be a way to significantly enhance signal without any further modification, providing better TOF approximation for low light level situations.