Translator Disclaimer
25 June 2018 Reconstruction of high-resolution early-photon tomography based on the first derivative of temporal point spread function
Author Affiliations +
For fluorescence molecular tomography, higher spatial resolution can be achieved using minimally scattered early photons. Conventional reconstruction methods of early photon tomography (EPT) are based on the integral of temporal point spread function (TPSF), which may lead to poor image quality due to systematic noise and time mismatch between the TPSF data and forward model. The derivative of the rising portion of TPSF is proposed to be used in EPT to increase the performance of reconstruction because the derivative is less sensitive to noise and time mismatch than the integral. A method based on projected Tikhonov regularization with the reconstructed result of steepest descent algorithm as a priori information is developed. Using the derivative of TPSF, the method can achieve high spatial resolution in phantom experiments and is capable of separating targets with an edge–edge distance of 1.5 mm.

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

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.2 In addition, it is reported that shorter time gate could lead to higher spatial resolution according to simulations.2 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:5

Eq. (1)

where c is the speed of light, r is the coordinate, and t is the time. ϕ(r,t) denotes the photon density, S(r,t) denotes the source of light, and D(r)=1/3(μs+μa), where μs is the reduced scattering coefficient and μa is the absorption coefficient.

Inverse problem of i-EPT can be formulated as

Eq. (2)

where Tgate is the time gate used. A(t)=[aij(t)] is the forward matrix, m is the number of discretized voxels, and n is the number of measurements. aij(t)=Wi(rj,t), where i denotes the i’th measurement, rj denotes the j’th discretized coordinate, and t is the time. In addition, Wi(r,t)=ϕim(r,t)*ϕie(r,t)*L(r,t) is convolution of time, where ϕim(r,t) denotes the sensitivity map of the detector of the i’th measurement, ϕie(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)=T0T0+Tgatefi(t)dt is derived from the i’th measurement corresponding to 0TgateA(t)dt, in which fi(t) denotes the TPSF of the i’th measurement and T0 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 T0 because of the limitation of temporal resolution of the detector. The time mismatch between the wrong T0 and the true T0 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:

Eq. (3)

where X denotes the fluorescence distribution. The i’th row of AS is given as AS(i)=dWi(r,t)/dt|t=Tgatei, where Tgatei is the time gate for the i’th measurement. YS(i)=dfi(T0+Tgatei)/dt is the slope of TPSF of the i’th measurement at the time t=T0+Tgatei corresponding to AS(i). In this letter, Tgatei=TpeakiT0, where Tpeaki is the time when the slope of the TPSF of the i’th measurement reaches its peak value and T0 is the starting point as mentioned before. To be noticed, each measurement has its own Tpeaki but the same starting point T0. Therefore, Tgatei varies with the measurements. A time mismatch ΔTmis will cause a ΔYS(i)=dfi(T0+ΔTmis+Tgatei)/dtdfi(T0+Tgatei)/dt in YS(i). ΔYS(i)=ΔTmisd2fi(T0+Tgatei)/dt2+o(ΔTmis) after Taylor’s expansion. When T0+Tgatei=Tpeak, d2fi(T0+Tgatei)/dt2=0, so ΔYS(i)=o(ΔTmis) is smaller. Therefore, the time gate Tgatei=TpeakiT0 is chosen.

Moreover, for numerical calculation, AS(i)=[Wi(r,Tgatei+Δt)Wi(r,TgateiΔt)]/2Δt and Ys(i)=k=1w[fi(Tgatei+kΔt)fi(TgateikΔt)]/w(w+1)Δt, where Δt=12.2  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×3×4  cm3 and discretization of 0.6×0.6×2  mm3. 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 l2-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:

  • Step 1: Steepest descent part,

    where Xk denotes the fluorescence distribution after k iterations, dk is the gradient, lk is the analytical step size, YS is the slope of TPSFs, and AS is the matrix derived from the forward matrix, as described previously. Iteration stops when k=K.K is the maximum number of iterations.

  • Step 2: Projected Tikhonov regularization,6

    XTKp+1=max{XTKp+[STS+trace(STS)αI]1ST(YSSXTKp),0},S={dA(i,j)dt},s.t.  XK(j)>0,
    where XTKp denotes the fluorescence distribution after p iterations in the region of interest where XK>0, and XK is the final reconstruction result of the steepest descent part. S denotes the matrix corresponding to XTKp. α is the regularization parameter. Furthermore, the initial value of the projected Tikhonov part is also set as the result of the steepest descent part.

Simulations and phantom experiments were based on the system built by our laboratory previously.7 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.7 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.8 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 μs=1  mm1 and absorption coefficient μa=0.003  mm1, 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Δt(Δt=12.2  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.9 For reconstruction, 800 iterations were used in the steepest descent algorithm and 3000 iterations were used in Tikhonov regularization, α=2×104. 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Δt, s-EPT is more tolerant to noise and time mismatch and can achieve better image quality and higher spatial resolution.

Fig. 1

Reconstruction results of (a) i-EPT and (b) s-EPT in simulations with 2% Gaussian noise and 5Δt mismatch added. (c) PCCs between data without noise and data affected with time mismatch from 5Δt to 5Δt and 2% Gaussian noise added to TPSFs.


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Δt to 5Δ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Δ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  μ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 α. 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 α was set to be a geometric progression with the first term as 102, common ratio as 0.667, and the last term no <104. The maximum k=8000 and the minimum α=1×104 were chosen empirically because the increase in k hardly affects the reconstruction result when k>8000 and an α<1×104 would probably lead to nonconvergent result. For each combination of k and α, 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 α. 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.

Fig. 2

PCCs-cross talk map of reconstruction results of s-EPT, i-EPT and all-time-gate FMT (a) EED=1.5  mm, (b) EED=3  mm, and (c) EED=6  mm.


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 EED=6  mm, k=200 and α=1×103; when EED=3  mm, k=1100 and α=5×104; when EED=1.5  mm, k=1700 and α=3×104. For i-EPT, when EED=6  mm, k=400 and α=103; when EED=3  mm, k=3200 and α=3×104; when EED=1.5  mm, k=6000 and α=3×104. For all-time-gate FMT, when EED=6  mm, k=200 and α=1×103; when EED=3  mm, k=1200 and α=5×104; when EED=1.5  mm, k=1800 and α=4×104. 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.

EED (mm)s-EPTi-EPTAll-time-gate

Fig. 3

Reconstruction results of (a1)–(a3) s-EPT, (b1)–(b3) i-EPT, and (c1)–(c3) all-time-gate FMT. (d1) Profiles along the white dashed lines in (a1)–(c1). (d2) Profiles along the white dashed lines in (a2)–(c2). (d3) Profiles along the white dashed lines in (a3)–(c3).


Although all three methods perform almost the same for EED=6  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.9 s-EPT is also found to perform better. Spatial resolution of separating targets with EED=1.5  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 background10,11 and this method might achieve good in vivo performance.


The authors have no relevant financial interests in this letter and no potential conflicts of interest to disclose.


This work was supported by the National Natural Science Foundation of China under Grant Nos. 81471665 and 81561168023.



V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods, 7 (8), 603 –614 (2010). 1548-7091 Google Scholar


F. Leblond et al., “Early-photon fluorescence tomography: spatial resolution improvements and noise stability considerations,” J. Opt. Soc. Am. A, 26 (6), 1444 –1457 (2009). JOAOD6 0740-3232 Google Scholar


Y. Mu and M. Niedre, “A fast SPAD-based small animal imager for early-photon diffuse optical tomography,” in 36th Annual Int. Conf. of the IEEE, Engineering in Medicine and Biology Society (EMBC), 2833 –2836 (2014). Google Scholar


L. Zhang et al., “Early-photon guided reconstruction method for time-domain fluorescence lifetime tomography,” Chin. Opt. Lett., 14 (7), 071702 (2016). CJOEE3 1671-7694 Google Scholar


B. Zhang et al., “Early-photon fluorescence tomography of a heterogeneous mouse model with the telegraph equation,” Appl. Opt., 50 (28), 5397 –5407 (2011). APOPAI 0003-6935 Google Scholar


A. Neumaier, “Solving ill-conditioned and singular linear systems: a tutorial on regularization,” SIAM. Rev., 40 (3), 636 –666 (1998). SIREAD 0036-1445 Google Scholar


C. Cai et al., “Direct reconstruction method for time-domain fluorescence molecular lifetime tomography,” Opt. Lett., 40 (17), 4038 –4041 (2015). OPLEDP 0146-9592 Google Scholar


C. Cai et al., “Nonlinear greedy sparsity-constrained algorithm for direct reconstruction of fluorescence molecular lifetime tomography,” Biomed. Opt. Exp., 7 (4), 1210 –1226 (2016). BOEICL 2156-7085 Google Scholar


J. Riley et al., “Choice of data types in time resolved fluorescence enhanced diffuse optical tomography,” Med. Phys., 34 (12), 4890 –4900 (2007). MPHYA6 0094-2405 Google Scholar


L. Sinha et al., “Getting more early photons with less background: detection rate and signal-to-background improvements in enhanced early photon imaging,” Proc. SPIE, 10486 104860E (2018). PSISDG 0277-786X Google Scholar


L. Sinha et al., “Enhanced detection of early photons in time-domain optical imaging by running in the” dead-time” regime,” Opt. Lett., 41 (14), 3225 –3228 (2016). OPLEDP 0146-9592 Google Scholar
© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2018/$25.00 © 2018 SPIE
Jiaju Cheng, Chuangjian Cai, and Jianwen Luo "Reconstruction of high-resolution early-photon tomography based on the first derivative of temporal point spread function," Journal of Biomedical Optics 23(6), 060503 (25 June 2018).
Received: 27 February 2018; Accepted: 5 June 2018; Published: 25 June 2018

Back to Top