Computationally-efficient linear scheme for overlap time-gating spatial frequency domain diffuse optical tomography using an analytical diffusion model

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.


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][2][3][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][6][7][8][9][10][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][6][7][8][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][14][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][22][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.

Forward problem
The forward problem of TD-SFD-DOT is solved by applying the TD-SFD-pDE.The fluencerate 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: Φ(r, t, f x ) = Φ 0 (r, t, f x )exp(iω x x + φ), with ω x = 2πf x the angular frequency and spatial phase φ of the SFD modulation.Φ 0 (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 Φ 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: where c is the speed of light in the medium; µ a and µ ′ s are the absorption and reduced scattering coefficients of the medium, respectively; D = 1/3/(µ a + µ ′ s ) is the diffusion coefficient.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 • k/t) and k = (4Dc) −1 is shown as follows: where z 1 = z − z 0 , z 2 = z + z 0 + 2z b , and z 0 = (µ a + µ ′ s ) −1 .z is the vertical distance to the medium surface, z b = D/A is the extrapolation distance with A = (1 − R eff )/(2 + 2R eff ) and R eff represents an effective reflection coefficient [26,27].The flux amplitude Γ 0 at the surface of the medium can be calculated using the formula Γ 0 (t, f x ) = AΦ 0 (z = 0, t, f x ).

Inverse problem
To solve the inverse problem of TD-SFD-DOT, the discretized matrix equation for solving µ a and µ ′ s perturbations is obtained from TD-SFD-pDE, as follows: where δΓ 0 is the vector containing the wide-field TR flux amplitude perturbations, δµ a and δµ ′ s are vectors that describe the absorption and reduced scattering perturbations at each inner node of the tissue domain, J a and J s represent the absorption and the component of scattering Jacobian matrices, respectively.To obtain the analytical form of J a and J s , the Green's function in semi-infinite geometry derived from TD-SFD-pDE is used: where ρ is the horizontal distance from detection location r d to the source location r = [ρ, z] and z is the vertical distance from r to the medium surface.
The superscript (Γ) denotes that Eq. ( 4) represents the flux on the medium surface.Equation (5) and Eq. ( 6) provide the analytical expressions for J a and J s elements at time t, respectively: where V r is the volume of the voxel at r, and the gradient operator in Eq. ( 6) is simplified to d/dz as Φ 0 varies only along the z-axis.The abbreviations used in the above expressions are displayed as: 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 τ 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 Φ 0 , G 0 (Γ) 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 h1 and h2 : where , and erf represents the error function.The variation of Jacobian matrices J a (r d , r z | ρ=0 , t, f x ) and J s (r d , r z | ρ=0 , t, f x ) are shown in Fig. 1, which represent the matrix elements of absorption and scattering at a given r d along the z-direction when ρ = 0, where the spatial frequency of the matrices is 0.02 mm −1 , µ a and µ ′ s are 0.008 mm −1 and 0.8 mm −1 , respectively.
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 Φ 0 and G 0 (Γ) that tend to 0 near the surface and the mirror sources are δΓ 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: , 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π, 4/3π 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 ) • 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 J a ⊗ IRF and J s ⊗ IRF, where ⊗ denotes the convolution operator.

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 ∆T b is divided into multiple overlap time bins of width ∆T g and a temporal center interval ∆T.A time bin centered at t n is obtained by averaging all TR samples over the time-gating range [t n − ∆T g /2, t n + ∆T g /2], giving the time-gating form of I(t) [19]: where M = ∆T g /∆T b represents the number of samples in the time-gating range.Equation (3) in its overlap time-gating form is shown as: where δΓ ∆T g 0 , J ∆T g a , and J ∆T g s represent the overlap time-gating form of δΓ 0 , J a ⊗ IRF and J s ⊗ IRF, respectively.
To analyze the temporal distribution of the absorption and scattering perturbation sensitivities, we calculate Jacobian matrices with intensity correction using expressions Ĵa (r d , r, and J ∆T g s (r d , r, t n , f x ) represent the values of Jacobian matrices at each node r in the medium for any sampling point 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 Γ ∆T g 0 (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 Ĵa and Ĵ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 r d and r.The vertical axes represent the values of Ĵa and Ĵ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 Ĵa with time, but it remains relatively stable, suggesting a more consistent distribution of absorption sensitivity.Compared with Ĵa , Ĵ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.
Since Jacobian matrices for overlap time-gating TD-SFD-DOT have four dimensions: spatial frequencies f x , wide-field sampling points r d , time gates T g , and inner nodes of the medium r, yielding the total number of the matrix elements N f x × N r d × N T g × N 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 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 ⟨r d ⟩ 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.
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, δµ (1)   a and δµ ′ (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 δµ (2)  a was discarded due to the selected early time gates are not sensitive to absorption and the resulting output: δµ a = δµ (1)  a and δµ ′ s = δµ ′(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.

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

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.
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 . U i (i = 1, 2) represent the optical properties of the Target 1 and 2, and U ′ 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 Target 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.

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 , 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 ′ 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 ′ .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.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 ′ 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.   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 ′ and k ′ 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.

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.
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 × N r d × N 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.

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.
Disclosures.The authors declare that there are no conflicts of interest related to this article.

Fig. 1 .
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.

Fig. 2 .
Fig. 2. Spatial variation of absorption and scattering perturbation sensitivities at different temporal centers and spatial frequencies.(a) The variation of Ĵa along the z-axis (ρ = 0 mm).(b) The horizontal variation of Ĵ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 Ĵs .

Fig. 3 .
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.

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.

Fig. 6 .
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.

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

Fig. 8 .
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 (
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 ′ presented inTable 2 demonstrate that the individual reduced scattering coefficient

Fig. 9 .
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.