Open Access
Add to BibSonomy
Abstract
Time-domain (TD) spatial frequency domain (SFD) diffuse optical tomography (DOT) potentially enables laminar tomography of both the absorption and scattering coefficients. Its full time-resolved-data scheme is expected to enhance performances of the image reconstruction but poses heavy computational costs and also susceptible signal-to-noise ratio (SNR) limits, as compared to the featured-data one. We herein propose a computationally-efficient linear scheme of TD-SFD-DOT, where an analytical solution to the TD phasor diffusion equation for semi-infinite geometry is derived and used to formulate the Jacobian matrices with regard to overlap time-gating data of the time-resolved measurement for improved SNR and reduced redundancy. For better contrasting the absorption and scattering and widely adapted to practically-available resources, we develop an algebraic-reconstruction-technique-based two-step linear inversion procedure with support of a balanced memory-speed strategy and multi-core parallel computation. Both simulations and phantom experiments are performed to validate the effectiveness of the proposed TD-SFD-DOT method and show an achieved tomographic reconstruction at a relative depth resolution of ∼4 mm.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Spatial frequency domain (SFD) imaging is a non-contact technique that provides rapid wide-field imaging of tissue optical properties using a sinusoidal modulated illumination pattern projected onto the tissue surface. Quantitative 2D maps of the absorption and scattering contrasts at different spatial frequencies [1–4] are used to achieve depth resolution, which is not a tomographic method to obtain the 3D distribution of optical properties. SFD diffuse optical tomography (SFD-DOT) enables tomographic imaging without dense source-detector measurements. Currently, SFD-DOT systems with continuous wave (CW) illumination have been developed, aiming at solving the absorption coefficient of the turbid medium [5–11]. CW-SFD-DOT reconstruction scheme is obtained either by matching the measurement images in the SFD and the DOT scheme in the spatial domain with the Fourier method [5–9], or by direct derivation with the phasor diffusion equation (pDE) [10,11]. However, the absence of temporal information measurement in the CW-SFD-DOT system causes difficulty in separating absorption and scattering contrasts.
To solve this problem, time-resolved (TR) measurement has been incorporated into the framework of SFD-DOT, as the temporal profiles offer the richest information about the tissue and have the potential to separate absorption and scattering perturbations [12,13]. SFD-DOT system in the time domain (TD) applies a wide-field temporally ultrashort pulsed source with spatial modulation, and the re-emitted photons from the tissue typically persist for a few nanoseconds in practice [13–15]. The propagation of the broad SFD patterns in the TD has been demonstrated, and the use of time-gating data allows the decoupling of absorption and scattering contrasts in the turbid medium [15,16]. Fourier transform is used to derive the Green’s function for time-dependent SFD-DE in infinite space, indicating that the absorption coefficient of the medium is spatial-frequency-dependent [15]. Several time-gating detection SFD-DOT systems have been reported. The temporal point spread function (TPSF) of the medium is detected with an ultrafast gating charge-coupled device camera [13,14], and the application of narrower gate widths is demonstrated to reduce the crosstalk of absorption and scattering perturbations [13]. For the purpose of reconstructing the optical properties of the tissue, Monte Carlo (MC) simulation is used to produce time-gating Jacobian matrices in accordance with illumination patterns [13,14,17], which provides high accuracy but is computationally intensive when simulating a large number of photons. Time-gating TD-SFD-DOT schemes described above employ a limited number of time bins, resulting in larger gate intervals and a sparse temporal sampling that does not fully exploit the information in the TR data. To fully utilize temporal information, it is necessary to employ measurements with higher temporal resolution in imaging systems such as time correlated single photon counting (TCSPC)-modules, which possess resolution of picoseconds. In addition, time-gating Jacobian matrices in the SFD have four dimensions containing information on structured light fields, wide-field sampling points, time gates, and inner nodes of the medium, which result in high computational resource consumption when solving the inverse problem. To reduce computational complexity, an approach based on generalized pulse spectrum technique has been extended to TD-SFD-DOT to characterizing the temporal profile by Laplace transform [12], saving the computation time by applying featured data, but the dimensionality reducing operation of temporal data leads to an increased level of effort required in decoupling absorption and scattering coefficients.
Full TR data is capable of achieving better absorption and scattering contrasts, although the datatype still requires overcoming the low signal-to-noise ratio (SNR) of the TCSPC-based measurement data. Dividing the TR data into multiple time bins can be an effective method [18], but it may result in a reduction in temporal resolution as the number of time bins decreases for better SNR, which is contrary to the purpose of setting narrower gate widths to reduce the crosstalk of optical properties. To balance the improvement of robustness to noise with the preservation of temporal information, overlap time-gating approach is a viable strategy that averages wide-field temporal data within the overlapping time-gating ranges [19]. Because of the temporal dispersion of impulse response function (IRF) in the system, it is necessary to match the measured TR data and the forward model predictions through the convolution of system's IRF with TR data predicted by the forward model [19,20], or using IRF to solve for the true photon time-of-flight distributions (TOFDs) [21–23]. Furthermore, the utilization of full TR information in the overlap time-gating TD-SFD-DOT makes it essential to fully adapt available computational resources to handle the large Jacobian matrices efficiently.
In this paper, we present a computationally-efficient overlap time-gating scheme of TD-SFD-DOT for tomographically solving the absorption and scattering contrasts in the turbid medium, which fully utilizes temporal information in the TR data, with the goal of obtaining wide-field optical properties in the shallow part of biological tissue. We adopt the analytical solution of TD-SFD-pDE in semi-infinite geometry as forward model for TD-SFD-DOT, and derive Jacobian matrices analytically for the reconstruction with overlap time-gating data that enhances the robustness of the scheme to noise while reducing the loss of temporal information in full TR data. To better reconstructing absorption and scattering perturbations, we develop a two-step linear inversion procedure utilizing full TR data. Due to the computational complexity of the inverse problem caused by the large size of Jacobian matrices, we employ the algebraic reconstruction technique (ART) [24,25] in the inversion procedure to avoid manipulating the entire Jacobian matrices. Meanwhile, we use multi-core parallel computation and adopt a strategy of balancing memory and speed to reconstruct perturbations, thereby maximizing the use of practically-available resources and improving solving efficiency. To validate the presented TD-SFD-DOT scheme for laminar tomography of optical properties, simulations and phantom experiments are performed. We further investigate the depth resolution of the proposed method and discuss the effect of the overlap time gate setting strategy on reconstructions.
2. Methods
2.1 Forward problem
The forward problem of TD-SFD-DOT is solved by applying the TD-SFD-pDE. The fluence-rate in the turbid medium stimulated by temporally pulsed spatial pattern illumination with a spatial frequency of ${f_x}$ in the x direction can be described in phasor form as: $\Phi (\mathbf{r},t,{f_x}) = {\Phi _0}(\mathbf{r},t,{f_x})\textrm{exp}(i{\omega _x}x + \varphi )$, with ${\omega _x} = 2\pi {f_x}$ the angular frequency and spatial phase $\varphi $ of the SFD modulation. ${\Phi _0}(\mathbf{r},t,{f_x})$ is the temporal and spatial distribution of fluence-rate amplitude stimulated by the AC component of the TD-SFD source. Inserting $\Phi $ into TD-DE and simplifying the equation gives the TD-SFD-pDE as Eq. (1), which describes the amplitude distribution of the fluence-rate along the z-axis in a semi-infinite homogeneous medium:
By applying the extrapolated boundary condition, the time-dependent fluence-rate amplitude distribution in the semi-infinite homogeneous medium stimulated by TD-SFD source with abbreviations $f(l) = \exp ( - {{l \cdot k} / t})$ and $k = {(4Dc)^{ - 1}}$ is shown as follows:
2.2 Inverse problem
To solve the inverse problem of TD-SFD-DOT, the discretized matrix equation for solving ${\mu _a}$ and ${\mu ^{\prime}_s}$ perturbations is obtained from TD-SFD-pDE, as follows:
Equation (5) and Eq. (6) are derived with extrapolated boundary condition, which introduces mirror photon sources that lead to inaccuracies in the fluence-rate description near the time 0, and therefore affect the elements of Jacobian matrices. To prevent this issue, we change the integration limits of each element in the Jacobian matrices on the time variable $\tau $ from $(0,t)$ to $({t_0},t - {t_0})$, where ${t_0}$ is positive and much smaller than t. With this process, the calculation of ${\Phi _0}$, ${G_0}^{(\Gamma )}$ and their derivatives with respect to z as time tends to 0 is avoided. Abbreviations ${h_1}$ and ${h_2}$ in Eq. (5) and Eq. (6) are replaced by ${\tilde{h}_1}$ and ${\tilde{h}_2}$:
Fig. 1. Jacobian matrices vary along the z-axis at a given detection point with horizontal distance of 0 mm and a spatial frequency of 0.02 mm-1. (a) and (b) show the variation of ${J_a}$ and ${J_s}$, and (c) shows the positive value of ${J_s}$. The black solid lines represent the elements with the integral interval $(0,t)$, while the red dashed lines represent those with the integral interval $({t_0},t - {t_0})$ at a ${t_0}$ of 0.1 ps.
Figure 1(a) and Fig. 1(b) show the variation of ${J_a}$ and ${J_s}$ along the z-axis, and Fig. 1(c) shows the positive values of ${J_s}$ along the z-axis, indicating that the abrupt numerical changes caused by the time variables of ${\Phi _0}$ and ${G_0}^{(\Gamma )}$ that tend to 0 near the surface and the mirror sources are smoothed. Furthermore, as shown in Fig. 1(b) and Fig. 1(c), the absolute values of the negative values of ${J_s}$ in the near-surface region are much larger than the positive values of ${J_s}$, leading to the concentration of the solved perturbations of ${\mu _a}$ and ${\mu ^{\prime}_s}$ at shallow depths. All negative values in the scattering Jacobian matrices are set to 0 empirically to balance the values in the matrices, enabling the resolution of deeper tomographic images.
$\delta {\mathbf{\Gamma }_0}$ in Eq. (3) is obtained using the Born approximation. To calibrate the absolute source intensity of the system and match the forward model prediction, a reference measurement for a homogeneous medium with known optical properties is performed. The measured data is calibrated by taking the ratio of the SFD-modulated amplitude distributions of the target and reference medium, expressed as: ${I_r}(t,{f_x}) = {I_{tar}}(t,{f_x})/{I_{ref}}(t,{f_x})$, where ${I_{tar}}$ represents the amplitude distribution on the surface of the target medium stimulated by TD-SFD illumination and demodulated through three measurements with pattern offsets $0$, $2/3\pi $, $4/3\pi $ rad, respectively. ${I_{ref}}$ represents the amplitude distribution on the reference medium surface, following the same measurement and demodulation method. Therefore, the calibrated TR amplitude distribution of the target medium is ${I_r}(t,{f_x}) \cdot {I_{pred}}(t,{f_x})$, where ${I_{pred}}$ represents the forward model prediction. Meanwhile, Eq. (3) requires system calibration by convolving time-dependent Jacobian matrices with IRF, as ${\mathbf{J}_a} \otimes IRF$ and ${\mathbf{J}_s} \otimes IRF$, where ${\otimes} $ denotes the convolution operator.
2.3 Reconstruction method
To enhance the robustness of the reconstruction while maximizing the utilization of temporal information, an overlap time-gating approach described in [19] is employed in wide-field TR data processing. The measured TR data $I(t)$ with a sampling interval $\Delta {T_b}$ is divided into multiple overlap time bins of width $\Delta {T_g}$ and a temporal center interval $\Delta T$. A time bin centered at ${t_n}$ is obtained by averaging all TR samples over the time-gating range $[{t_n} - \Delta {T_g}/2,{t_n} + \Delta {T_g}/2]$, giving the time-gating form of $I(t)$ [19]:
To analyze the temporal distribution of the absorption and scattering perturbation sensitivities, we calculate Jacobian matrices with intensity correction using expressions ${\hat{\boldsymbol{J}}_a}({{\mathbf{r}_d},\mathbf{r},{t_n},{f_x}} )= {{\boldsymbol{J}_a^{\Delta {T_g}}({{\mathbf{r}_d},\mathbf{r},{t_n},{f_x}} )} / {\Gamma _0^{\Delta {T_g}}({{t_n},{f_x}} )}}$ and ${\hat{\boldsymbol{J}}_s}({{\mathbf{r}_d},\mathbf{r},{t_n},{f_x}} )= {{\boldsymbol{J}_s^{\Delta {T_g}}({{\mathbf{r}_d},\mathbf{r},{t_n},{f_x}} )} / {\Gamma _0^{\Delta {T_g}}({{t_n},{f_x}} )}}$. $\boldsymbol{J}_a^{\Delta {T_g}}({{\mathbf{r}_d},\mathbf{r},{t_n},{f_x}} )$ and $\boldsymbol{J}_s^{\Delta {T_g}}({{\mathbf{r}_d},\mathbf{r},{t_n},{f_x}} )$ represent the values of Jacobian matrices at each node $\mathbf{r}$ in the medium for any sampling point ${\mathbf{r}_d}$ at z = 0 mm, respectively. The absorption and reduced scattering coefficients for matrix calculation are 0.008 mm-1 and 0.8 mm-1 respectively. After IRF-calibration, Jacobian matrices are converted into their overlap time-gating form with ${t_n}$ as the temporal center and a gate width of 400 ps. To ensure a meaningful comparison, all Jacobian matrices are convolved with the same IRF, and the matrices are divided by IRF-convolved $\Gamma _0^{\Delta {T_g}}({{t_n},{f_x}} )$ to ensure that the sensitivity of absorption and scattering are not influenced by the time-varying intensity of the TPSF.
The spatial variations of ${\hat{\boldsymbol{J}}_a}$ and ${\hat{\boldsymbol{J}}_s}$ at different temporal centers and spatial frequencies are displayed in Fig. 2. The horizontal axes labels ρ and z, respectively, denote the horizontal and vertical distances of ${\mathbf{r}_d}$ and $\mathbf{r}$. The vertical axes represent the values of ${\hat{\boldsymbol{J}}_a}$ and ${\hat{\boldsymbol{J}}_s}$. The sensitivity of absorption shown in Fig. 2(a) and Fig. 2(b) appears to be higher at the late time range of measurements, as indicated by the increase in the value of ${\hat{\boldsymbol{J}}_a}$ with time, but it remains relatively stable, suggesting a more consistent distribution of absorption sensitivity. Compared with ${\hat{\boldsymbol{J}}_a}$, ${\hat{\boldsymbol{J}}_s}$ decreases rapidly with time, as plotted in Fig. 2(c) and Fig. 2(d), indicating that scattering perturbation is more sensitive in the early time range, and therefore the solution for scattering perturbation can be obtained by primarily using TR data within the early time range. In the late time range of the TCSPC-based measurements, which are more sensitive to absorption perturbation, the low SNR caused by the reduction of photon count leads to an increase in the difficulty of reconstructing perturbations by applying data only in the late time range, which indicates the need to use full TR information in absorption perturbation reconstruction.
Fig. 2. Spatial variation of absorption and scattering perturbation sensitivities at different temporal centers and spatial frequencies. (a) The variation of ${\hat{\boldsymbol{J}}_a}$ along the z-axis (ρ = 0 mm). (b) The horizontal variation of ${\hat{\boldsymbol{J}}_a}$ with respect to ρ (z = 1.5 mm), where the negative signs of ρ indicate that the direction of spatial variation is opposite to that of positive values. (c) and (d) display the corresponding variation of ${\hat{\boldsymbol{J}}_s}$.
Since Jacobian matrices for overlap time-gating TD-SFD-DOT have four dimensions: spatial frequencies ${f_x}$, wide-field sampling points ${\mathbf{r}_d}$, time gates ${T_g}$, and inner nodes of the medium $\mathbf{r}$, yielding the total number of the matrix elements ${N_{{f_x}}} \times {N_{{\mathbf{r}_d}}} \times {N_{{T_g}}} \times {N_\mathbf{r}}$, it is impractical to manipulate the entire absorption and scattering Jacobian matrices simultaneously due to memory limitations. To overcome this problem, we apply ART to reconstruct optical contrasts, which allows row-by-row computations without loading the entire matrices.
Figure 3 shows the flowchart that illustrates the computationally efficient process we use to solve for perturbations based on ART. For a specific measurement system with a fixed reference phantom, IRF is certain and the analytical-calculated Jacobian matrices are not updated during the ART inversion process. Therefore, Jacobian matrices are computed and saved before optical contrast reconstruction to reduce the amount of computation. TD-Jacobian matrices for each sampling point ${\mathbf{r}_d}$ are obtained by computing Eq. (5) and Eq. (6). IRF-calibration is then performed in the dimension t of the matrices before they are converted into their overlap time-gating form and saved. This process is executed in parallel on a multi-core computer to increase the efficiency of calculation. For the ART inversion, wide-field sampling points on the medium surface are divided into n groups randomly, with the number of groups adapted to the practical computing resources available to balance memory and speed for efficient computation. Jacobian matrices for each sampling point set ${\left\langle {{\mathbf{r}_d}} \right\rangle _i}$ are loaded into memory and removed after finishing the corresponding part of ART inversion, then the next set of Jacobian matrices is loaded for reconstruction until all sampling points are traversed and an ART inversion is completed.
Fig. 3. Flowchart of the Jacobian matrix computation and single ART-based inversion process for overlap time-gating TD-SFD-DOT. Square brackets indicate the dimensions of the corresponding matrices.
Furthermore, to improve the separation of absorption and scattering and achieve better contrasts, a two-step inversion procedure utilizing overlap time-gating data in different time ranges is presented by considering the difference in the temporal distribution characteristics of absorption and scattering sensitivities. As shown in Fig. 4, the first linear inversion in this procedure employs time-gating data for the entire range in the TD with initial perturbations of 0. Its solutions, $\delta \mu _a^{(1)}$ and $\delta \mu ^{\prime {(1)}}_s$, are set as initial values for the second inversion, which uses early time gates to better separate the absorption and scattering perturbations. The early time gates are defined as the gates with temporal centers preceding the peak time of TPSFs in the reference measurement. The solved absorption perturbation $\delta \mu _a^{(2)}$ was discarded due to the selected early time gates are not sensitive to absorption and the resulting output: $\delta {\mu _a} = \delta \mu _a^{(1)}$ and $\delta {\mu ^{\prime}_s} = \delta \mu ^{\prime {(2)}}_s$. For both steps, ART linear inversions are performed according to the flowchart shown in Fig. 3. A 24-core computer (Intel Xeon CPU E5-2650 v4, 2.20 GHz) is utilized to solve for optical perturbations using absorption and scattering Jacobian matrices each with 2 × 2500 × 301 × 25000 elements, which are then converted to their overlap time-gating form resulting in 2 × 2500 × 34 × 25000 elements. It takes ∼56 min to build Jacobian matrices and ∼27 min to execute an entire two-step inversion procedure for solving optical contrasts, with the first step taking ∼22 min and the second employing early overlap time-gate dataset and taking ∼5 min.
3. Results and discussion
3.1 Experimental setup
To validate the effectiveness of the proposed overlap time-gating TD-SFD-DOT method, simulations and phantom experiments are performed. The experiment system and the phantom structure are shown in Fig. 5.
Fig. 5. TR-SFD-DOT measurement system and target phantom setting. (a) Experimental of TR-SFD measurement; (b) Structure of the target phantom; (c) Photo of the target phantom illuminated at a spatial frequency of 0.02 mm-1, with three-phase projected patterns.
As shown in Fig. 5(a), a 780 nm pulsed laser is modulated by a digital micro-mirror device (DMD) to produce wide-field sinusoidal illumination patterns with a spatial frequency of 0.02 mm-1. The reflected wide-field temporal data is obtained using single-pixel imaging method, which utilizes a detection DMD for spatial sampling and a single-photon avalanche diode (SPAD)-based TCSPC module for TR measurement with a sampling interval of 12.5 ps. The illumination field is 60 × 60 (mm) and the detection range is 50 × 50 (mm). A 40 mm thick homogeneous liquid phantom with an absorption coefficient of 0.008 mm-1 and a reduced scattering coefficient of 0.8 mm-1 constructed using Intralipid (20% fat emulsion injection, Fresenius Kabi, Jiangsu, China) and India Ink (PH1714, Phygenge Biotechnology, Fuzhou, China) is used for the reference measurement. For the target phantom given in Fig. 5(b), two 12 mm diameter tubes, labeled Target 1 and Target 2, filled with Intralipid-Ink mixture that provides optical contrasts are placed in the homogeneous phantom. The absorption and scattering contrasts between Target 1 and the phantom background are both 3:1, while for the Target 2, the corresponding contrasts are both 2:1. Additionally, the photo of the target phantom with SFD illumination is provided in Fig. 5(c), and the picture also displays illumination patterns. The IRF of the system is measured by replacing the phantom with an opaque, reflective slab.
3.2 Experimental validation
The tomographically reconstructed perturbations in the phantom experiments are presented in Fig. 6(a), where the DC and AC component at spatial frequency 0.02 mm-1 are combined, the width of the overlap time gates is set to 400 ps with a temporal interval of 100 ps, and the proposed two-step inversion procedure is used for reconstruction. Figure 6(b) displays the corresponding perturbation reconstructions performed using the ART-based typical one-step inversion procedure. For comparison with the phantom experiments, photon TOFDs on the surface of the target phantom in Fig. 5(b) are generated using MC simulation platform [28]. Figure 6(c) shows the optical property reconstructions obtained by applying the time-gating data, which has the same setup as the phantom experiment and is converted from the simulated TOFDs that convolved with the system’s IRF. Figure 6(d) displays the perturbation reconstructions using the overlap time gates which are obtained from generated TOFDs without IRF convolution, and have a gate interval of 20 ps and a gate width of 50 ps. First eight overlap gates are used for the reconstruction, and the inverse problem is solved using typical one-step inversion procedure based on ART because of the monotonic decreasing and the high SNR of the simulated TOFDs.
Fig. 6. Reconstructed images of absorption and reduced scattering coefficients. (a) shows the phantom experiment reconstructions in the x-y plane at z = 2.75 mm and in the x-z plane at y = 25 mm and the line profiles along the x-axis at z = 2.75 mm obtained by the proposed two-step inversion procedure, while (b) displays the corresponding solutions with the typical one-step inversion procedure. (c) and (d) present reconstruction images using simulated TR data with and without IRF convolution, respectively. The actual positions of the heterogeneities are indicated by dashed lines.
The cross-sections and x-line profiles of the reconstructed absorption and reduced scattering coefficients depicted in Fig. 6 show the effectiveness of the proposed TD-SFD-DOT scheme for separating absorption and scattering contrasts. To facilitate a quantitative comparison, two evaluation factors are defined as follows: relative error of the maximum value ${Q_m} = |{U_i} - \textrm{max(}{U^{\prime}_i})|/{U_i} \times 100\%$ and contrast ratio ${Q_r} = \textrm{max(}{U^{\prime}_1})/\textrm{max(}{U^{\prime}_2})$. ${U_i}(i = 1,2)$ represent the optical properties of the Target 1 and 2, and ${U^{\prime}_i}(i = 1,2)$ denote the peak of the corresponding reconstructions in the x-profile of at a given depth. Accordingly, the theoretical value of ${Q_r}$ is 1.5. The quantitative evaluation of the x-profiles at z = 2.75 mm for each case in Fig. 6, is presented in Table 1.
Table 1. Quantitative evaluation of reconstructions shown in Fig. 6
Compared to the reconstructions obtained with the typical inversion procedure based on ART, the utilization of our proposed two-step procedure provides a better quantification of the optical property perturbations, as evidenced by the overall improvement in ${Q_m}$ and ${Q_r}$, especially for scattering. In addition, it is observed that the absorption reconstructions obtained with the proposed method are more consistent with the actual shape of the inclusions. A comparison of phantom experimental and MC-simulated cases, both of which are affected by the IRF, reveals that the overall difference in the solved values of perturbations results in ${Q_m}$ fluctuations. Furthermore, the factors ${Q_r}$ for the IRF-convolved simulation data demonstrate an increase in the reconstructed absorption and scattering contrasts. Considering the measurement noise and errors in the experimentally collected TR data, the above phenomenon is reasonable. In the case of MC-simulation without IRF-convolution, the utilization of the generated TOFDs in the early time range results in a loss of absorption information, which may in turn lead to a depth offset and large ${Q_m}$. However, the absorption- and scattering-related ${Q_r}$ for this case are most closely aligned with their expectations, while other cases influenced by IRF exhibit contrast reduction, indicating the impact of temporal information dispersion caused by IRF. In addition, vertical cross-sections shown in Fig. 6 indicate that the depth of tomographically reconstructing perturbations are limited, which might be due to the decay of the SFD illumination along the z-axis, even though a low-spatial-frequency modulated illumination is used and the inverse problem is solved by combining the DC component. Therefore, only laminar perturbations within a depth of 5 mm are solved.
3.3 Depth resolution validation
To discuss the depth resolution of the proposed TD-SFD-DOT scheme, a series of simulation experiments are performed. The simulation TR data are generated using the MC simulation platform [28]. The target medium contains a cylinder heterogeneity with a radius of 1 mm, placed in the center of the medium near the surface and oriented in the same direction as the cylinder shown in Fig. 5(b). The heterogeneity has absorption and reduced scattering coefficients twice those of the background, and the background absorption and reduced scattering coefficients are 0.008 mm-1 and 0.8 mm-1, respectively. The depth d from the axis of the heterogeneity to the medium surface is set to 1.5, 2, 2.5 and 3 mm. To quantify the depth of the reconstructed heterogeneities, the centroid depth, defined as $d^{\prime} = (\sum\nolimits_i^Zu_i^{\prime} z_i)/(\sum\nolimits_i^Z u_i^{\prime} )$, is computed from the normalized z-profile of the reconstruction at x = 25 mm and y = 25 mm, where Z is the number of nodes, ${u^{\prime}_i}$ is the normalized reconstruction value of the i-th node and ${z_i}$ denotes the depth of the node. The centroid depth distributions are linearly fitted, and the relative depth resolution of each case is compared based on the fitted slope $k^{\prime}$. The tomographic images shown in Fig. 7(a) and Fig. 7(b) are reconstructed using simulated temporal data without IRF convolution, while images shown in Fig. 7(c) and Fig. 7(d) are reconstructed using the temporal data that convolved with IRF. To better display the reconstructions of the central heterogeneity, the edges of the profiles on the x-z plane are discarded. The centroid depths and fitted slopes for the above cases are presented in Table 2. The solutions indicate that our proposed method can achieve simultaneous reconstruction of the absorption and scattering contrasts at a depth resolution of ∼4 mm without IRF convolution, despite some depth offset. However, the IRF-calibration in the DOT inverse problem makes it difficult to decouple absorption and scattering perturbations due to the dispersion of temporal information caused by the IRF-convolution, which leads to a lack of depth resolution.
Fig. 7. Verification of depth resolution for the simultaneous reconstruction of absorption and scattering perturbations. The reconstructed (a) absorption and (b) reduced scattering profiles using simulated TR data without IRF convolution on the x-z plane at y = 25 mm, and the corresponding normalized line profiles along z-axis at x = 25 mm and y = 25 mm. (c) and (d) present the absorption and reduced scattering reconstruction images using simulated TR data with IRF convolution. The actual positions of the heterogeneities are indicated by dashed lines.
Table 2. Depth resolution evaluation for each case using MC-simulated TR data
Furthermore, the depth resolution of the proposed scheme is discussed when the absorption and scattering perturbations are solved separately with IRF-convolved simulation data. The simulation employs the same structure and the background optical properties. When solving for the absorption perturbation only, the absorption contrast of the heterogeneity is 2:1 with the background, and the reduced scattering coefficient is equal to the background. Similarly, when solving for the scattering perturbation only, the heterogeneity has a scattering contrast of 2:1 with the background, and its absorption coefficient is equivalent to the background. The simulated TR data convolved with IRF is converted into overlap time gates with an interval of 100 ps and a width of 400 ps, and the inversion is implemented using typical ART.
Figure 8(a) displays the results of utilizing full time-gating data to solve the absorption perturbation, indicating that the proposed scheme can achieve a relative depth resolution of absorption contrast within ∼4 mm. Compared to solutions without IRF convolution in Fig. 7(a), results in Fig. 8(a) show enhanced accuracy at the gap between the heterogeneity and the medium surface, accompanied by higher $k^{\prime}$ as indicated in Table 2. This may be attributed to different time gate selection strategies, which suggests that the use of full temporal information through the convolution between IRF and the simulated TR data leads to a better tomographic reconstruction for absorption. The process of the solutions depicted in Fig. 7(a) utilizes the generated TOFDs in the early time range, resulting in a loss of absorption information. In contrast, the process of reconstructing results presented in Fig. 8(a) employs full temporal information through the convolution between IRF and the simulated TR data, leading to a better tomographic reconstruction for absorption.
Fig. 8. Verification of depth resolution for separately solving absorption and scattering perturbations using IRF-convolved simulated TR data. (a) Absorption profiles that reconstructed with full temporal data on the x-z plane at y = 25 mm, and the corresponding normalized line profiles along z-axis at x = 25 mm and y = 25 mm. (b) and (c) show the reduced scattering profiles that reconstructed using simulated temporal data in the full and early time range respectively. The actual positions of the heterogeneities are indicated by dashed lines.
Figure 8(b) and Fig. 8(c) show the images of scattering perturbation using full and early time-gating data, respectively, indicating that the introduction of IRF leads to a degradation of depth resolution in the results obtained by solving only for the scattering perturbation. The fitted slopes $k^{\prime}$ presented in Table 2 demonstrate that the individual reduced scattering coefficient reconstructions, both with full and early time gates, are capable of providing relative depth resolution. It should be noted that the presence of an abnormal value in the z-profile of the full-TR-data-based reconstruction for the target medium containing the heterogeneity at a depth of 3 mm, results in a shift of the centroid depth. The calculated $d^{\prime}$ and $k^{\prime}$ without the abnormal value are supplemented in Table 2 and labelled with parentheses. The reconstructions using early TR data are observed to be more stable with fewer artifacts and overall higher relative depth resolution performance, suggesting that utilizing early time gates to solve scattering perturbations is reasonable in the proposed two-step inversion procedure. In comparison to the case of solving with true TOFDs, the degradation of depth resolution might be due to the dispersion of temporal information about the scattering contrast carried by photons emitted in the early time range caused by IRF. Meanwhile, the convolution performed to calibrate IRF in the reconstruction scheme also disperses the temporal information in the Jacobian matrices, making it difficult to solve the depth of the scattering perturbation in the inverse problem.
In general, although the proposed TD-SFD-DOT inversion scheme improves the separation of optical contrasts, it still requires further enhancement as the dispersion of TPSFs induced by IRF results in a reduction of depth resolution, indicating a necessity for strengthening the capacity to extract information from TPSFs. Therefore, it is necessary to apply other IRF-calibration methods, or to recover the true TOFDs from TR measurements. The degradation of depth resolution might also be related to the IRF pulse width, which is beyond the scope of this paper. Furthermore, the incorporation of additional parameters related to optical perturbation sensitivities, such as spatial frequency, into the inversion procedure has the potential to enhance the depth resolution while concurrently reducing the crosstalk of absorption and scattering.
3.4 Time gate setup discussion
The effect of varying time gate setting strategies on the solved perturbations by applying the TCSPC-based measurement data in the phantom experiments is presented in Fig. 9, and the corresponding evaluation factors ${Q_m}$ and ${Q_r}$ for the x-profile with depth z = 2.75 mm are listed in Table 3. Figure 9(a)-(d) show tomographic reconstructions solved with various gate widths and a fixed temporal interval of 100 ps. The absorption reconstruction using non-overlapping time gates with both gate width and interval of 100 ps shows less accuracy in the description of inclusion shapes, resulting in increased ${Q_m}$ and the factor ${Q_r}$ with the most significant offset. As the gate width increases, the time range of gates begins to overlap, which facilitates an increase in signal quality by averaging more time samples. Meanwhile, as the center interval of time gates is fixed, the temporal resolution of the overlap time-gating data remains constant, preventing the degradation of temporal information and resulting in perturbation solutions with improved image quality and quantification. However, increasing the gate width further degrades the quality of the reconstructions. When the gate width is increased to 600 ps, a general decline in reconstruction performance is observed in both qualitative and quantitative factors, which might be related to different SNRs at each sampling point of the TCSPC-based measured TPSFs. A moderate gate width will reduce the noise by integrating over the temporal data. If there is too much overlap between the gates, the photon counting data in the low-SNR range will be combined with the data in the high-SNR range and converted to time gates, leading to a loss of data quality in the high-SNR range. In addition, since the scattering contrasts are solved in two steps and with TR data in the early time range, their reconstructions are less affected by the time gate setup than the absorption contrasts, resulting in more stable reconstructions.
Fig. 9. Reconstruction profiles on the x-y plane at z = 2.75 mm and the x-z plane at y = 25 mm with the corresponding line profiles along x-axis at z = 2.75 mm. The reconstructed (a)-(b) absorption and (c)-(d) reduced scattering coefficient using different gate widths with the same temporal interval of 100 ps, and the reconstructed (e)-(f) absorption and (g)-(h) reduced scattering coefficient using different intervals with the same gate width of 400 ps. The actual positions of the heterogeneities are indicated by dashed lines.
Table 3. Performance comparison of reconstructions presented in Fig. 9.
Figure 9(e)-(h) present absorption and reduced scattering perturbation reconstructions with different time gate intervals and a fixed gate width 400 ps. As the temporal interval increases, the contrast and quantification of the solved optical properties decrease due to the redundancy of temporal information, indicating that a smaller gate interval enhances the reconstruction quality. However, in the process of solving the inverse problem, a decrease in gate interval leads to a larger matrix equation to be inversed. Adding each additional time gate results in an increase of ${N_{{f_x}}} \times {N_{{\mathbf{r}_d}}} \times {N_\mathbf{r}}$ elements in the each Jacobian matrix. Therefore, depending on the solution accuracy and computational efficiency requirements, it is necessary to set the overlap time gate appropriately.
4. Conclusion
In this paper, we present a computationally-efficient overlap time-gating scheme of TD-SFD-DOT to solve the tomographic image of absorption and scattering perturbations simultaneously. We adopt the analytical solution of TD-SFD-pDE in semi-infinite geometry, from which the forward model is derived and analytically formulated the Jacobian matrices in the overlap time-gating form to enhance the SNR of the TR data while reducing redundancy and fully utilizing the temporal information. Additionally, we propose a two-step linear inversion procedure based on ART to improve the ability to separate absorption and scattering contrasts. Multi-core parallelism and a balanced memory-speed strategy are employed for efficient computation by utilizing available computing resources in the inversion process. Simulations and phantom experiments are performed to validate the proposed scheme, demonstrating its capability to solve tomographic images of absorption and scattering within a depth of ∼4 mm. The depth resolution of the proposed TD-SFD-DOT scheme with IRF-calibration is discussed, with the conclusion that the IRF-caused dispersion to TPSFs leads to a degradation of depth resolution. This highlights the necessity of further enhancing the capacity of separating optical contrasts. This paper also discusses the impact of the time gate setting strategy on the quality of the reconstructed images, showing that selecting a moderate gate width and reducing the gate interval can enhance the stability of the solution.
The proposed scheme in this paper will be applied to wide-field shallow biological tissue tomography, with further expansion to wide-field fluorescence molecular tomography in the TD, where we expect to have a broader range of applications. Future work and enhancement will focus on applying more precise models for structured light propagation in the TD-SFD-DOT scheme and employing more effective IRF-calibration methods with the aim of preventing the dispersal of temporal information. In addition, spatial frequency and other potential parameters will be incorporated into the inversion procedure to achieve enhanced absorption and scattering contrasts.
Funding
National Natural Science Foundation of China (62075156, 62205239, 81971656, 82071971).
Disclosures
The authors declare that there are no conflicts of interest related to this article.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. David J. Cuccia, F. Bevilacqua, A. J. Durkin, et al., “Modulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain,” Opt. Lett. 30(11), 1354–1356 (2005). [CrossRef]
2. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, et al., “Quantitation and mapping of tissue optical properties using modulated imaging,” J. Biomed. Opt. 14(2), 024012 (2009). [CrossRef]
3. Rolf B. Saager, David J. Cuccia, and A. J. Durkin, “Determination of optical properties of turbid media spanning visible and near-infrared regimes via spatially modulated quantitative spectroscopy,” J. Biomed. Opt. 15(1), 017012 (2010). [CrossRef]
4. T. A. Erickson, A. Mazhar, D. Cuccia, et al., “Lookup-table method for imaging optical properties with structured illumination beyond the diffusion theory regime,” J. Biomed. Opt. 15(3), 036013 (2010). [CrossRef]
5. S. D. Konecky, A. Mazhar, D. Cuccia, et al., “Quantitative optical tomography of sub-surface heterogeneities using spatially modulated structured light,” Opt. Express 17(17), 14780–14790 (2009). [CrossRef]
6. V. Lukic, V. A. Markel, and J. C. Schotland, “Optical tomography with structured illumination,” Opt. Lett. 34(7), 983–985 (2009). [CrossRef]
7. A. Bassi, C. D’Andrea, G. Valentini, et al., “Detection of inhomogeneities in diffusive media using spatially modulated light,” Opt. Lett. 34(14), 2156–2158 (2009). [CrossRef]
8. C. D’Andrea, N. Ducros, A. Bassi, et al., “Fast 3D optical reconstruction in turbid media using spatially modulated light,” Biomed. Opt. Express 1(2), 471–481 (2010). [CrossRef]
9. S. Bélanger, M. Abran, X. Intes, et al., “Real-time diffuse optical tomography based on structured illumination,” J. Biomed. Opt. 15(1), 016006 (2010). [CrossRef]
10. T.-X. Li, Z.-P. Qin, X. Hou, et al., “Multi-wavelength spatial frequency domain diffuse optical tomography using single-pixel imaging based on lock-in photon counting,” Opt. Express 27(16), 23138–23156 (2019). [CrossRef]
11. M. Dan, T.-X. Li, X. Hou, et al., “A lock-in photon-counting based single pixel imaging toward real-time multi-wavelength SFD-DOT,” Proc. SPIE 11234, 112341 L (2020). [CrossRef]
12. W.-X. Bai, Y.-H. Dong, Y.-R. Zhang, et al., “Wide-field illumination diffuse optical tomography within a framework of single-pixel time-domain spatial frequency domain imaging,” Opt. Express 32(4), 6104–6120 (2024). [CrossRef]
13. V. Venugopal, J. Chen, and X. Intes, “Development of an optical imaging platform for functional imaging of small animals using wide-field excitation,” Biomed. Opt. Express 1(1), 143–156 (2010). [CrossRef]
14. J. Chen, V. Venugopal, F. Lesage, et al., “Time-resolved diffuse optical tomography with patterned-light illumination and detection,” Opt. Lett. 35(13), 2121–2123 (2010). [CrossRef]
15. A. Bassi, C. D’Andrea, G. Valentini, et al., “Temporal propagation of spatial information in turbid media,” Opt. Lett. 33(23), 2836–2838 (2008). [CrossRef]
16. V. Venugopal, J. Chen, F. Lesage, et al., “Time resolved diffuse optical tomography using a digital light processor,” in Biomedical Optics and 3-D Imaging (Optica Publishing Group, 2010), p. BSuD52.
17. R.-Y. Yao, X. Intes, and Q.-Q. Fang, “Direct approach to compute Jacobians for diffuse optical tomography using perturbation Monte Carlo-based photon “replay,” Biomed. Opt. Express 9(10), 4588–4603 (2018). [CrossRef]
18. H.-J. Zhao, F. Gao, Y. Tanikawa, et al., “Time-resolved diffuse optical tomography and its application to in vitro and in vivo imaging,” J. Biomed. Opt. 12(6), 062107 (2007). [CrossRef]
19. J. Li, X. Yi, X. Wang, et al., “Overlap time-gating approach for improving time-domain diffuse fluorescence tomography based on the IRF-calibrated Born normalization,” Opt. Lett. 38(11), 1841–1843 (2013). [CrossRef]
20. W.-B. Wan, L.-L. Liu, Y.-H. Wang, et al., “Time-resolved early-photon scheme for high-resolution fluorescence molecular tomography with perturbation Monte Carlo modeling: Experimental validations using a multichannel TCSPC system,” Infrared Phys. Technol. 98, 323–333 (2019). [CrossRef]
21. M. Diop and K. St. Lawrence, “Deconvolution method for recovering the photon time-of-flight distribution from time-resolved measurements,” Opt. Lett. 37(12), 2358–2360 (2012). [CrossRef]
22. D. A. Smith, G. McKenzie, A. C. Jones, et al., “Analysis of time-correlated single photon counting data: a comparative evaluation of deterministic and probabilistic approaches,” Methods Appl. Fluoresc. 5(4), 042001 (2017). [CrossRef]
23. J. T. Giurleo and D. S. Talaga, “Global fitting without a global model: Regularization based on the continuity of the evolution of parameter distributions,” J. Chem. Phys. 128(11), 114114 (2008). [CrossRef]
24. R. Gordon, R. Bender, and G. T. Herman, “Algebraic Reconstruction Techniques (ART) for three-dimensional electron microscopy and X-ray photography,” J. Theoret. Biol. 29(3), 471–481 (1970). [CrossRef]
25. X. Intes, V. Ntziachristos, J. P. Culver, et al., “Projection access order in algebraic reconstruction technique for diffuse optical tomography,” Phys. Med. Biol. 47(1), N1–N10 (2002). [CrossRef]
26. R. C. Haskell, L. O. Svaasand, T.-T. Tsay, et al., “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11(10), 2727–2741 (1994). [CrossRef]
27. A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14(1), 246–254 (1997). [CrossRef]
28. L.-M. Yu, F. Nina-Paravecino, D. Kaeli, et al., “Scalable and massively parallel Monte Carlo photon transport simulations for heterogeneous computing platforms,” J. Biomed. Opt. 23(01), 010504 (2018). [CrossRef]
