Fluorescence molecular tomography (FMT) is a vital modality capable of detecting in vivo biological activities on the molecular level. Thus, FMT is popular in cancer study, drug development, and other biological areas.¹ But due to strong scattering in biological tissues, FMT inverse problem is extremely illposed, which results in poor image quality, mainly low spatial resolution. To improve the spatial resolution of FMT, several methods have been developed, such as applying anatomical a priori information, using minimally scattered (snake-like) early photons and so on. According to previous simulations, spatial resolution can be significantly improved when using early arriving photons instead of full-time-gate photons as data to reconstruct images as early photons experience fewer scattering events.²

Conventional reconstruction methods of early photon tomography (EPT) are generally based on the integral of temporal point spread function (TPSF), which represents the total amount of photons received by a detector before a certain time gate.² In addition, it is reported that shorter time gate could lead to higher spatial resolution according to simulations.² In phantom experiments, resolution of integral-based EPT (i-EPT) is typically limited at separating targets with edge–edge distance (EED) about 6 mm.^3,4 But in phantom experiments, i-EPT behaves worse than all-time-gate tomography that uses the integral of all photons received and is similar to steady-state FMT. It may be because of high-level noise in early photons part or time mismatch between TPSF and forward model due to the limitation of temporal resolution of the detector, which will be both mentioned later in this letter.

In forward model calculation and simulations, telegrapher equation was applied, shown as follows:⁵

where $c$ is the speed of light, $r$ is the coordinate, and $t$ is the time. $\mathrm{\varphi }(r,t)$ denotes the photon density, $S(r,t)$ denotes the source of light, and $D(r)=1/3({\mu }_s^{\prime }+{\mu }_a)$, where ${\mu }_s^{\prime }$ is the reduced scattering coefficient and ${\mu }_a$ is the absorption coefficient.

Inverse problem of i-EPT can be formulated as

where $T_{\mathit{\text{gate}}}$ is the time gate used. $A(t)=[a_{ij}(t)]$ is the forward matrix, $m$ is the number of discretized voxels, and $n$ is the number of measurements. $a_{ij}(t)=W_i(r_j,t)$, where $i$ denotes the $i$’th measurement, $r_j$ denotes the $j$’th discretized coordinate, and $t$ is the time. In addition, $W_i(r,t)={\mathrm{\varphi }}_i^m(r,t)*{\mathrm{\varphi }}_i^e(r,t)*L(r,t)$ is convolution of time, where ${\mathrm{\varphi }}_i^m(r,t)$ denotes the sensitivity map of the detector of the $i$’th measurement, ${\mathrm{\varphi }}_i^e(r,t)$ denotes the excitation map of the $i$’th measurement, calculated by Eq. (1). $L(r,t)$ is the lifetime map of fluorescence molecule, which is usually unknown and set as an estimated value. $X$ denotes the fluorescence distribution to be reconstructed. $Y(i)={\int }_{T_0}^{T_0+T_{\text{gate}}}{f}_i(t)\mathrm{d}t$ is derived from the $i$’th measurement corresponding to ${\int }_0^{T_{\text{gate}}}A(t)\mathrm{d}t$, in which $f_i(t)$ denotes the TPSF of the $i$’th measurement and $T_0$ is the time point when the laser excitation starts, i.e., the starting point. It should be noted that it is difficult to determine the exact $T_0$ because of the limitation of temporal resolution of the detector. The time mismatch between the wrong $T_0$ and the true $T_0$ will adversely affect the performance of i-EPT, as will be demonstrated later in this letter.

To improve the performance of EPT, reconstruction based on the first-order derivative (i.e., slope) of the rising portion of TPSF is proposed in this letter. The method is called as slope-EPT (s-EPT). Projected Tikhonov regularization is employed to solve the inverse problem. s-EPT needs to solve the following equation:

where $X$ denotes the fluorescence distribution. The $i$’th row of $A_S$ is given as $A_S(i)=d{W}_i(r,t)/dt{|}_{t=T_{\text{gate}}^i}$, where $T_{\text{gate}}^i$ is the time gate for the $i$’th measurement. $Y_S(i)=d{f}_i(T_0+T_{\text{gate}}^i)/dt$ is the slope of TPSF of the $i$’th measurement at the time $t=T_0+T_{\text{gate}}^i$ corresponding to $A_S(i)$. In this letter, $T_{\text{gate}}^i=T_{\text{peak}}^i-T_0$, where $T_{\text{peak}}^i$ is the time when the slope of the TPSF of the $i$’th measurement reaches its peak value and $T_0$ is the starting point as mentioned before. To be noticed, each measurement has its own $T_{\text{peak}}^i$ but the same starting point $T_0$. Therefore, $T_{\text{gate}}^i$ varies with the measurements. A time mismatch $\mathrm{\Delta }{T}_{\mathrm{mis}}$ will cause a $\mathrm{\Delta }{Y}_S(i)=d{f}_i(T_0+\mathrm{\Delta }{T}_{\mathrm{mis}}+T_{\text{gate}}^i)/dt-d{f}_i(T_0+T_{\text{gate}}^i)/dt$ in $Y_S(i)$. $\mathrm{\Delta }{Y}_S(i)=\mathrm{\Delta }{T}_{\mathrm{mis}}{d}^2{f}_i(T_0+T_{\text{gate}}^i)/d{t}^2+o(\mathrm{\Delta }{T}_{\mathrm{mis}})$ after Taylor’s expansion. When $T_0+T_{\text{gate}}^i=T_{\text{peak}}$, $d^2{f}_i(T_0+T_{\text{gate}}^i)/d{t}^2=0$, so $\mathrm{\Delta }{Y}_S(i)=o(\mathrm{\Delta }{T}_{\mathrm{mis}})$ is smaller. Therefore, the time gate $T_{\text{gate}}^i=T_{\text{peak}}^i-T_0$ is chosen.

Moreover, for numerical calculation, $A_S(i)=[W_i(r,T_{\text{gate}}^i+\mathrm{\Delta }t)-W_i(r,T_{\text{gate}}^i-\mathrm{\Delta }t)]/2\mathrm{\Delta }t$ and $Y_s(i)=\sum _{k=1}^w[f_i(T_{\text{gate}}^i+k\mathrm{\Delta }t)-f_i(T_{\text{gate}}^i-k\mathrm{\Delta }t)]/w(w+1)\mathrm{\Delta }t$, where $\mathrm{\Delta }t=12.2\mathrm{ps}$ is the temporal resolution of the photomultiplier tube (PMT) and $w$ is the length of a smoothing window for derivation and is empirically set as 20. Forward model was calculated within a cylinder with a diameter of 3 cm and height of 4 cm through the finite-element method. The result was interpolated into uniform mesh with a size of $3\times 3\times 4{\mathrm{cm}}^3$ and discretization of $0.6\times 0.6\times 2\mathrm{mm}$³. Discretization in the height direction was set as 2 mm mainly for calculation efficiency. During reconstruction, only a 0.8-cm-thick layer around the excitation plane was taken into consideration, also for calculation efficiency.

In this letter, a projected Tikhonov regularization algorithm combined with analytical-step-size steepest descent algorithm is proposed to solve the inverse problem given as Eq. (3). With analytical step size, steepest descent algorithm will be easily trapped in local optimum, whereas it will also be tolerant to noise. According to the stimulations, the projected Tikhonov regularization algorithm performs much better than the conventional Tikhonov regularization algorithm without projection, because the former is optimized to overcome the over-smoothing limitation of the $l_2$-norm regularization. With results from the steepest descent algorithm as a priori information, the projected Tikhonov regularization algorithm could achieve both high stability and high spatial resolution.

The algorithm can be described as follows:

Simulations and phantom experiments were based on the system built by our laboratory previously.⁷ The system is based on PMT. During measurements, the phantom was rotated for 360 deg with increments of 20 deg to obtain 18 projections. For each projection, 11 detectors were uniformly distributed within a detecting field of view (FOV) of 220 deg.⁷ A group of filters including an achromatic doublet (AC254-030-B, Thorlabs, Newton, New Jersey) and two bandpass fluorescence filters with a center wavelength of 840 nm (ff01-840/12-25, Semrock, Rochester, New York; XBPA840, Asahi Spectra, Torrance, California) were added to eliminate excitation light.⁸ In phantom experiments, femtosecond laser generator (Spectra-Physics, Newport Corporation, Canada) was set to work at 780-nm wavelength (80 MHz, 100-fs pulse-width).

In simulations, two homogeneous cylindrical targets with diameters of 4 mm were buried inside a 3-cm-diameter cylinder with reduced scattering coefficient ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and absorption coefficient ${\mu }_a=0.003{\mathrm{mm}}^{-1}$, which are the same as the coefficients of 1% Intralipid. Two targets were placed symmetrically with an EED of 3 mm. Time gate for i-EPT was set as the average of time gates for s-EPT. Moreover, to test the performance of both methods, a $-5\mathrm{\Delta }t$($\mathrm{\Delta }t=12.2\mathrm{ps}$ as mentioned above) time mismatch and 2% Gaussian noise were added to TPSFs from the measurements. Gaussian noise was added because in PMT system, noise could be considered mainly Gaussian.⁹ For reconstruction, 800 iterations were used in the steepest descent algorithm and 3000 iterations were used in Tikhonov regularization, $\alpha =2\times {10}^{-4}$. These parameters were empirically chosen. It could be seen from Figs. 1(a) and 1(b) that with 2% Gaussian noise and time mismatch of $-5\mathrm{\Delta }t$, s-EPT is more tolerant to noise and time mismatch and can achieve better image quality and higher spatial resolution.

To illustrate the reason, the Pearson correlation coefficient (PCC) between the original data without noise and the data affected by noise and time mismatch was calculated. In simulations, 2% Gaussian noise was added to the TPSFs generated and time mismatch varied from $-5\mathrm{\Delta }t$ to $5\mathrm{\Delta }t$. For each time mismatch, five data samples were generated. As shown in Fig. 1(c), the PCCs of s-EPT are close to 1, whereas those of i-EPT vary from 0.998 to 0.999. As can be seen from Fig. 1, with $-5\mathrm{\Delta }t$ time mismatch, although the PCC of i-EPT is only slightly lower than that of s-EPT with difference of $<0.002$, the reconstruction result of s-EPT is significantly better than that of i-EPT. It is mainly because of the ill condition of FMT inverse problem and the higher PCC of s-EPT explains why it performs better than conventional i-EPT in simulations. Moreover, as can be seen from the error bars (i.e., standard deviations), the data of s-EPT are much more stable than that of i-EPT and thus increases the robustness of reconstruction. To evaluate the possible spatial resolution of s-EPT, i-EPT and all-time-gate FMT, more stimulations with smaller EEDs were carried out. In this letter, we defined spatial resolution as the smallest EED with which the two targets could be separated. Two targets were considered separable when the cross talk of the two targets was lower than half maximum of the targets and the area of each reconstructed target was no smaller than half of the true area. According to these stimulations, the spatial resolution of s-EPT is about 0.4 mm, whereas those of i-EPT and all-time-gate FMT are about 1.2 and 1.1 mm, respectively. However, to be noticed, phantom experiments would be more complicated because noises are not only from the system and time mismatch but also from error of the forward model and other uncontrolled factors.

To further evaluate the performance of s-EPT, three sets of phantom experiments were carried out. Two tubes filled with $10\mu \mathrm{M}$ Indocyanine green/dimethylsulphoxide (ICG/DMSO) were buried symmetrically inside a 3-cm-diameter cylinder filled with 1% Intralipid. The tubes have inner diameters of 4 mm and outer diameters of 5 mm. The EEDs of the two tubes were 6, 3, and 1.5 mm, respectively. It should be noticed that EED in this letter refers to the distance of the inner edges of the tubes. All experiments were reconstructed through s-EPT, i-EPT, and all-time-gate FMT.

For all the methods, 800 iterations were used in the steepest descent algorithm. In the Tikhonov regularization algorithm, reconstruction results would be influenced by both the number of iterations $k$ and regularization parameter $\alpha$. To test the performance of the methods with different parameters, $k$ was set to range from 100 to 8000 with an increment of 100, whereas $\alpha$ was set to be a geometric progression with the first term as ${10}^{-2}$, common ratio as 0.667, and the last term no $<{10}^{-4}$. The maximum $k=8000$ and the minimum $\alpha =1\times {10}^{-4}$ were chosen empirically because the increase in $k$ hardly affects the reconstruction result when $k>8000$ and an $\alpha <1\times {10}^{-4}$ would probably lead to nonconvergent result. For each combination of $k$ and $\alpha$, the normalized cross talk defined as the minimum value along the central line of the two targets, and the PCC between the reconstruction result and the true values is calculated. Figure 2 shows the PCCs-cross talk map of the reconstruction results of different methods at all combinations of $k$ and $\alpha$. It can be observed that in experiments with EEDs of 1.5 and 3 mm, s-EPT is better than the other two methods as the points of s-EPT are typically closer to point (1,0), which means higher resolution and better image quality. In experiment with EED of 6 mm, the three methods performed similarly and could achieve good image quality and separate the two targets completely.

According to Fig. 2, the optimum parameters could be found to achieve both high PCC and low cross talk between targets. For s-EPT, when $\mathrm{EED}=6\mathrm{mm}$, $k=200$ and $\alpha =1\times {10}^{-3}$; when $\mathrm{EED}=3\mathrm{mm}$, $k=1100$ and $\alpha =5\times {10}^{-4}$; when $\mathrm{EED}=1.5\mathrm{mm}$, $k=1700$ and $\alpha =3\times {10}^{-4}$. For i-EPT, when $\mathrm{EED}=6\mathrm{mm}$, $k=400$ and $\alpha ={10}^{-3}$; when $\mathrm{EED}=3\mathrm{mm}$, $k=3200$ and $\alpha =3\times {10}^{-4}$; when $\mathrm{EED}=1.5\mathrm{mm}$, $k=6000$ and $\alpha =3\times {10}^{-4}$. For all-time-gate FMT, when $\mathrm{EED}=6\mathrm{mm}$, $k=200$ and $\alpha =1\times {10}^{-3}$; when $\mathrm{EED}=3\mathrm{mm}$, $k=1200$ and $\alpha =5\times {10}^{-4}$; when $\mathrm{EED}=1.5\mathrm{mm}$, $k=1800$ and $\alpha =4\times {10}^{-4}$. Figure 3 shows the reconstruction results with the optimum parameters, normalized by individual maximum values. In this letter, the optimum parameters are defined as those that achieve the highest PCC with cross talk being $<0.1$. Moreover, as shown in Table 1, PCCs for the reconstruction results along the central line agree with the maps shown in Fig. 2.

Table 1

The PCCs of the phantom experiments.

Although all three methods perform almost the same for $\mathrm{EED}=6\mathrm{mm}$ as shown in Figs. 3(a1)–3(d1), it can be seen from Figs. 3(a2)–3(d2) and 3(a3)–3(d3) that results of s-EPT are closer to true values and have fewer artifacts than the other two methods. In addition, results of i-EPT and all-time-gate FMT are too sparse when achieving cross talk lower than 0.1 and the shapes of the reconstructed targets are distorted. This is the main reason why their PCCs are lower than those of s-EPT. Moreover, as shown in Figs. 3(a2) and 3(d3), the results of i-EPT have a false peak, which suggests that its cross talk is higher than those shown in Fig. 2. Thus, the spatial resolution of i-EPT is lower, as measured by the cross talk. It should be noted that all-time-gate FMT performs better than s-EPT in relative value accuracy, as shown in Figs. 3(a2) and 3(d3). Nevertheless, in general, s-EPT performs better than the other two methods.

In conclusion, with the first-order derivative (i.e., slope) of the rising portion of TPSF as data, a method based on projected Tikhonov regularization with a priori information from steepest descent algorithm is proposed for EPT reconstruction. The method is proved to be more tolerant to noise than conventional methods and can achieve better image quality and higher resolution than both i-EPT and all-time-gate FMT. In addition, s-EPT is compared with EPT based on the peak of the TPSF reported in a previous study.⁹ s-EPT is also found to perform better. Spatial resolution of separating targets with $\mathrm{EED}=1.5\mathrm{mm}$ was achieved and further experiments with closer targets need to be carried out to find out the spatial resolution limit in phantom experiments. In the future, automated selection of parameters and in-vivo experiments should be performed. According to preliminary stimulations with nonuniform medium, it is possible for s-EPT to keep good performance in in vivo experiments. Moreover, the detection of the early photon could be enhanced with less background^10,11 and this method might achieve good in vivo performance.