Non-stationary reconstruction for dynamic fluorescence molecular tomography with extended kalman filter.

Dynamic fluorescence molecular tomography (FMT) plays an important role in drug delivery research. However, the majority of current reconstruction methods focus on solving the stationary FMT problems. If the stationary reconstruction methods are applied to the time-varying fluorescence measurements, the reconstructed results may suffer from a high level of artifacts. In addition, based on the stationary methods, only one tomographic image can be obtained after scanning one circle projection data. As a result, the movement of fluorophore in imaged object may not be detected due to the relative long data acquisition time (typically >1 min). In this paper, we apply extended kalman filter (EKF) technique to solve the non-stationary fluorescence tomography problem. Especially, to improve the EKF reconstruction performance, the generalized inverse of kalman gain is calculated by a second-order iterative method. The numerical simulation, phantom, and in vivo experiments are performed to evaluate the performance of the method. The experimental results indicate that by using the proposed EKF-based second-order iterative (EKF-SOI) method, we cannot only clearly resolve the time-varying distributions of fluorophore within imaged object, but also greatly improve the reconstruction time resolution (~2.5 sec/frame) which makes it possible to detect the movement of fluorophore during the imaging processes.


Introduction
Dynamic fluorescence molecular tomography (FMT) is important for bio-medical research since it allows non-invasive, tomographic imaging of the dynamic bio-distributions of fluorescent bio-markers within small animals in vivo [1][2][3]. With the development of optical contrast agents [4,5], more applications can be expected in detecting drug efficacy, monitoring disease progression, and evaluating novel therapy, etc. However, challenges remain in the reconstruction of dynamic FMT. The main reason is the current reconstruction methods [6-9] mainly focus on solving the stationary FMT tomographic problems. They are performed based on the (implicit) assumption that the fluorophore is stationary that it does not change during FMT imaging processes. Generally, the imaging time is needed to collect a complete set of measurements, which will be used as input data of the conventional reconstruction methods. However, in some applications, e.g., the drug delivery research, the fluorophore's concentration and location may change significantly when collecting these projection data. If the conventional reconstruction methods are directly applied to these time-varying fluorescence measurements, the reconstructed results may suffer from a high level of artifacts [10,11]. In addition, based on these stationary methods, only one tomographic image can be obtained after scanning one circle projection data. Considering the rotation and the exposure time in in vivo experiments, the acquisition time is generally >1 min [3]. As a result, some fast changes of fluorophore in imaged object may not be detected during FMT imaging processes. To some extent, the above factors limit the widespread application of dynamic FMT in bio-medical research.
To address the problem, some theoretical and experimental solutions have been proposed [10][11][12][13][14]. Considering that during FMT imaging processes, the fluorophore may be assumed as a discrete time stochastic process, it is possible to solve the timing-varying fluorescence tomographic problems by the state-estimation approach [11][12][13]. As demonstrated by S. Prince and V. Kolehmainen et al. [11,12], using the state-estimation method combined with extend kalman filter (EKF) technique [15], they successfully imaged the kinetic behaviors of fluorophore within imaged object. Despite of these successful applications, the current reconstruction method based on kalman filter was demonstrated by the fiber-based systems of slab geometry. Compared with the widely-used full-angle FMT imaging system, the fiberbased systems have some limitations in spatial sampling and experimental simplicity [16]. In addition, in the current EKF methods, the kalman gain was calculated based on the direct inverse matrix of the forward model of FMT. However, in FMT, the forward model matrix is highly ill-posed due to the high degree of light scattering in biological tissue. Together, the kalman gain is a core step in kalman filter processes. As a result, the imaging performance of FMT obtained by kalman filter may be affected if the kalman gain is calculated by the direct inverse method.
To overcome these limitations, in this paper, we apply kalman filter to time-varying fluorescence measurements which are acquired from a hybrid free-space, full-angle FMT/xray computed tomography (XCT) imaging system. Especially, to improve the imaging performance of FMT obtained by kalman filter, the generalized inverse of kalman gain is calculated by a second-order iterative method. The simulation, phantoms, and in vivo experiments are performed to evaluate the performance of the proposed method. Briefly, in the numerical simulation, two cases are demonstrated. Firstly, based on a three-dimensional (3-D) digital mouse model [17], we simulate the time-varying bio-distributions of indocyanine green (ICG) in the mouse's kidneys. In the second case, we simulate the movement of the fluorophore in the 3-D digital mouse during FMT imaging processes. For the phantom study, one tube filled with ICG is placed in one glass cylinder, which is used to imitate the timing-varying fluorescence imaging processes. In the in vivo experiment, utilizing the hybrid imaging system, we further image the bio-distributions of ICG in a nude mouse liver after the intravenous injection. The experimental results indicate that by using the proposed EKF-based second-order iterative (EKF-SOI) method, we can clearly resolve the dynamic bio-distributions of fluorophore within imaged object. In addition, the reconstruction time resolution can be greatly reduced (~2.5 sec/frame) in contrast to the conventional stationary reconstruction method, which makes it possible to detect the movable fluorophore during the imaging processes.
The outline of this paper is as follows. In section 2, the methods are detailed. In section 3, we describe the experimental materials. In section 4, the simulation, phantom, and in vivo experimental results are shown. Finally, we discuss the results and draw conclusions in section 5.

FMT imaging model
For the continuous wave FMT with point excitation source, the fluorescence measurement, Φ , acquired from surface of the imaged object can be generated based on the first order Born approximation [18] ( , ) is the point excitation source. By discretization, Eq. (1) can be rewritten in following form where ( ) ⋅ ℘ is a forward model used to map the unknown fluorescence distribution to known measurements. In this paper, the diffusion equation is solved using COMSOL Multiphysics 3.3 (COMSOL, Inc., Burlington, MA, USA).

Reconstruction method based on EKF
For the conventional methods, in reconstruction processes, the fluorophore ρ does not change during the acquisition of the measurement vector Φ . However, in some applications, the fluorophore's concentration and location may change significantly during imaging processes, as mentioned above. If these conventional reconstruction stationary methods are used to solve the time-varying measurements, on one hand, the reconstructed results may suffer from a high level of artifacts [11,12]. On the other hand, only one tomographic image can be obtained after scanning one circle projection data. As a result, it is impossible to image the movement of fluorophore during the imaging processes.
To address the problem, the state space approach is used in this paper. Briefly, in FMT imaging, the fluorophore ρ can be assumed as a discrete time stochastic Markov process, as follows, where t ℑ is a state transition matrix; t n is a zero mean state noise process with a covariance matrix t Q ; t is a time steps (state index). In addition, the fluorescence measurements can be derived from t ρ using the following forward model, is a forward model (observation matrix) corresponding to t ρ ; t m is an observation noise with a covariance matrix t R . Equations (4) and (5) constitute the state-space representation of the time-varying FMT problem. In the linear case, the kalman filter, as a typical algorithm, can be used to solve the problem. However, in FMT, the forward model is generally a nonlinear function. For the nonlinear problem, the kalman filter is not applicable. To address this problem, the nonlinear Eq. (5) is firstly approximated to the following linear equation, ℘ is a Jacobian matrix. After that, the kalman filter method is used. This method is called the extended kalman filter. The full equations for the (extend) kalman filter is given as follows, State prediction : Variance prediction : Kalman gain : ( ) where t G is a kalman gain; t C is a estimation error covariance matrix for the kalman filter state vector; t Q is a covariance matrix associated with noise vector t n ; t R is a covariance matrix of the observation errors t m . Based on Eqs. (7)-(11), we find that the fluorophore t ρ is calculated from current and past data. In order to incorporate the whole data in the sequence, the fixed-interval kalman smoother is used, which is achieved by performing a backward run with the following equation ρ is the kalman smoother estimate and N is a total number of time points. The smoothing gain 1 t − ℵ is generated as follows,

Second-order iteration for kalmain gain
Let T = + A JCJ R , if the generalized inverse of A can be obtained, the kalmain gain can be calculated by Eq. (9). As mentioned above, in FMT, the forward model J is highly ill-posed. As a result, the imaging performance of FMT obtained by EKF may be affected, if the kalman gain is generated based on the direct inverse matrix of A . To more effectively calculate the generalized inverse of A , in this paper, a second-order iterative method [19,20] is used as follows, where 0 S is chosen as where T represents transposition operation. Referring to [19], the series will be convergent to the generalized inverse of A when k → ∞ . The detailed derivation can be found in [19].

Experimental setup
The numerical simulation, phantom, and in vivo experiments were performed based on a hybrid FMT/XCT imaging system [21]. As shown in Fig. 1, the imaged object was fixed on a rotating stage, which allowed rotation of the imaged object for full-angle projection acquisitions. Around the rotation stage, a full-angle, free-space FMT and micro-XCT system were constructed, respectively. The detailed information on the imaging system could be found in [21].

Numerical simulation
In the numerical simulation, the imaged object was generated from a digital mouse model [17], as shown in Fig. 2. The absorption ( a μ ) and scattering ( ' s μ ) coefficients from [22] were assigned to different organs of the mouse (liver, kidneys, stomach, spleen, and pancreas), respectively, to simulate photons propagation in heterogeneous tissues. The optical properties outside these organs were regarded as homogeneous. The intensity curve of ICG in mouse kidneys, which is acquired from [24]. Inset in (c) shows that ICG intensities for 24 fluorescence projections corresponding to 24 angles from circle 1 (acquired during 0 min to 1 min). Different colors correspond to actual boundaries of different organs (red: livers; green: kidneys; stomach: cyan; pancreas: magenta; spleen: blue; bones: yellow; white: surface).
For FMT imaging, the collimated light beam [see red points in Fig. 2(a)] was modeled as an isotropic point source, located at one mean free path of photon transportation beneath the surface, at the heights of 0.9 cm and 1.7 cm, respectively. In the imaging processes, o 360 full view imaging was performed with 24 projections in a 15 o step.

Case 1: simulated the metabolic processes of ICG in kidneys of the mouse
With the imaging model, we firstly simulated the time-varying metabolism processes of ICG in the mouse's kidneys. Briefly, according to Fig. 2(c), different ICG intensities were set to the kidneys at every time points during measurement data collection process. For example, when collecting the first circle of fluorescence projection data (acquired during 0-60s), 24 different ICG intensities were set to the kidneys at each measured time point ( 1 2 4 t t − ). After that, fluorescence projection images from different angles were generated from Eq. (1).
Here, ICG intensity curve was acquired from [24]. Based on our in vivo experiments [3,14], the interval of adjacent projection was set to 2.5 sec.

Case 2: simulated the movement of fluorophore in the mouse
In the case 2, we simulated the movement of fluorophore in the mouse during the imaging time ( 1 2 4 t t − ) which is needed to collect a complete set of data (e.g., 24 projection images in the case) used as the input data of the conventional reconstruction methods. In detail, as shown in Fig. 3, in 1 6 t t − , the fluorophore (see the black cylinder in Fig. 3) was assumed to locate in the position 1. From 7 t to 14 t , the fluorophore moved to the position 2. Finally, the fluorophore changed to the position 3 between 15 t to 24 t . Here, a cylinders (2 mm in diameter and 4 mm in height) were used as fluorescence targets. In the simulation processes, the optical properties and the imaging parameters are the same with those used in the case 1. In the imaging processes, we needed to collect a complete set of data (24 projection images), which will be used as the input data of the conventional reconstruction methods. Briefly, In 1 6 t t − , the fluorophore is assumed to locate the position 1. In 7 1 4 t t − , the fluorophore moves to the position 2. From 15 t to 24 t , the fluorophore changes to the position 3.

Phantom experiment
The phantom experiment was performed to further validate the performance of the proposed method. As shown in Fig. 4(a), the phantom was made of a glass cylinder (outer diameter of 3.0 cm) filled with a mixture of intralipid, water, and ink, with an absorption coefficient of  Measurement data were obtained from the hybrid FMT/micro-XCT imaging system. In FMT imaging process, 24 fluorescence images were collected at every o 15 . Subsequently, 72 white light images were acquired in o 5 steps to form a 3-D geometry of the imaged phantom. Finally, a steel anchor point, which could be imaged in both two imaging systems, was used to provide the relative height information for registration.
Referring to [23], in FMT imaging processes, the fluorescence measurement Φ is approximately linearly with the fluorescence yield ρ . So, for simplification, in the phantom experiment, during the acquisition of a circle of data, the time-varying measurements are generated based on the following equation where k Φ is the acquired fluorescent photon density of the th k projection angle in the phantom experiment; ˆk Φ is the synthetics fluorescent photon density of the th k projection angle; k ξ is the ICG concentration at the th k time point (angle), which is acquired from the concentration curve in Fig. 4(b). After obtaining the full-angle synthetics measurements (24 angles in the case), these data are assembled to form the time-varying fluorescence measurement sequence, which are used to as the input data of the reconstruction.

In vivo experiment
In the in vivo imaging, an eight-week-old nude mouse, injected with 100 L μ of 50 / g mL μ ICG via tail vein, was suspended on the rotation stage and then continuously rotated 120 circles to monitor the metabolism processes of ICG in the liver. The acquisition time for collecting one circle of projection data was about 1 min. For each circle, 24 fluorescence images were acquired in 15 o steps. To correct for optical inhomogeneities, 24 excitation light images were subsequently acquired.
After optical imaging, a hepatobiliary contrast agent (Fenestra LC) was slowly intravenously injected to enhance the anatomical structure of livers in XCT imaging. 60 minutes after the injection of the contrast agent, the mouse was scanned using the micro-XCT to obtain the anatomical structures of organs. Finally, similar to the phantom experiments, a steel anchor point was used to provide the registration information. Throughout experiments, anesthesia was maintained with an isoflurane veterinary vaporizer.

Simulation study
4.1.1. Reconstructed the metabolic processes of ICG through the kidneys of the mouse Figure 5 shows the reconstructed tomographic images, obtained by the conventional stationary reconstruction method and the proposed reconstruction method based EKF, respectively. Here, the 24 projection data from the first circle were used as the input data of the reconstruction methods. When collecting one frame data, the fluorophore's concentrations vary significantly. Figures 5(a) Figure 6 shows the reconstructed FMT image sequence at all time points ( 1 2 4 t t − , with a time interval 2.5s) during the data acquisition time for collecting the first circle of projection data (0-60s). Furthermore, Fig. 6(b) shows the 3-D rendering of the reconstructed FMT image sequence. Here, the entire tomographic image sequence (24 frame) is obtained by the proposed EKF-SOI method.
The results indicate that it is difficult to resolve the fluorescence distributions in the kidneys based on the conventional stationary method [see Figs. 5(a) and 5(d)]. In contrast, when extend kalman filter is used in reconstruction processes, we can clearly resolve the ICG bio-distributions in the kidneys, as shown in Figs. 5(b), 5(c), 5(e), and 5(f). Furthermore, it can be seen that in this case, the reconstruction results from the proposed EKF-DIM method are similar with those obtained by the EKF-SOI method. On the other hand, we can also find that by using the conventional reconstruction method, only one tomographic image [see Figs. 5(a) and 5(d)] can be obtained based on the acquired 24 projection data from the first circle. For comparison, based on the same data set, 24 tomographic images ( 1 2 4 t t − ) corresponding to different angle can be obtained by the proposed EKF-SOI reconstruction method (see Fig.  6). That is, based on the EKF-SOI method, we can obtain the entire tomographic images at all time points ( 1 2 4 t t − , with a time interval 2.5s) during FMT imaging processes. In other words, by using the proposed EKF-SOI method, the reconstruction time resolution of FMT is reduced from 1 min (the stationary reconstruction method) to ~2.5 sec, which will be helpful for observing the fast biological activities.  Figure 7 compares the tomographic images obtained by the stationary, EKF-SOI, and EKF-DIM methods, respectively. Especially, during the acquisition time (0-60s) that is needed to collect a complete set of measurements for traditional reconstruction method, the fluorophore moved three times within the mouse. The detailed information can be found in Section 3.2.2. Figure 7(a) shows the reconstruction result obtained by the stationary method. Figures 7(b)-7(e) and Figs. 7(f)-7(i) show the reconstruction result obtained by the EKF-SOI and EKF-DIM method, respectively. In the reconstruction processes, the initialization values of the EKF were set as 0 ρ 0 = , 0 0.025 = C I, 0 0.01 = Q I, and 0 0.0001 = R I and t ℑ was chosen as identity matrix. In the EKF-SOI reconstruction, the number of the second-order iteration was set as 20. These parameters were the same as those used in case 1 (Section 4.1.1).

Reconstructed the movement of fluorophore in the mouse
The results indicate that based on the acquired 24 projection data, only one tomographic image was generated by the stationary reconstruction method. As a result, it is impossible to resolve the different positions of fluorophore during imaging processes [see Fig. 7(a)]. In contrast, when the proposed EKF-SOI method is used in reconstruction processes, based on the same projection data, the entire 24 tomographic images at all time points (0-60 s, with a time interval 2.5 s) can be generated [see Figs In addition, similar to the case 1, the results also indicate that there is no significant difference between the EKF-SOI and EKF-DIM results (see the 2nd and 3rd rows of Fig. 7). Further, Fig. 8 shows the 3-D rendering of the reconstructed FMT image sequence, which is obtained by the proposed EKF-SOI method. Based on these reconstructed 3-D images ( 2.5s,15s, 35s, and 60s t = ), we can resolve the moveable fluorophore. For example, Fig. 8(b) indicates that the fluorophore locates in the position 1. Figure 8(c) indicates that the fluorophore moves to the position 2. Similarly, Fig. 8(d) indicates that the fluorophore further moves to the position 3. Fig. 8. The 3-D rendering of the reconstructed images at different time points (2.5s, 15s, 35s, and 60s). These reconstructed images are obtained by the proposed EKF-SOI reconstruction method based on one circle of projection data. All images are displayed in the same range. Figure 9 depicts the reconstruction images from the phantom experiment. Figure 9(a) shows the reconstruction results obtained by the stationary method, where one circle of projection data was used as the input data of the reconstruction. Figures 9(b) and 9(c) show the reconstruction results obtained by the EKF-SOI and EKF-DIM method, respectively, which were imaged at 17.5 s (corresponding to the 7th projection angle).

Phantom experiment
For the phantom reconstruction, the volume of interest was a 3.0 cm 3.0 cm 3.3 cm × × 3-D region and sampled to 31 31 17 × × voxels. 9,680 voxels inside the imaged object and 16,181 source-detector pairs were used for reconstruction. Similar to the numerical simulation, for the conventional stationary methods, the reconstruction was performed by ART method. In the EKF reconstruction, the parameters of the EKF were 0 ρ 0 = , 0 0.025 = C I, 0 0.01 = Q I, and 0 0.0001 = R I, respectively. For the EKF-SOI reconstruction, the second-order iteration was terminated after 20 iterations. Note that the time-varying measurements from 24 angles (one circle) were used as the input data of the reconstruction. Fig. 9. The reconstruction results from the phantom obtained by the conventional stationary and EKF methods, where the time-varying measurements from 24 angles were used as the input data of the reconstruction. (a) The reconstruction results obtained by the conventional method based on circle of projection data (24 projection data). (b) The reconstruction result obtained by the proposed EKF-SOI method, which is imaged at imaged at 17.5s (corresponding to the 7th projection angle). (c) The reconstruction result obtained by the EKF-DIM method, imaged at 17.5s. The green curves depict the phantom boundary obtained by back-projecting the 72 white light images [25], which are used to validate the accuracy of registration. Figure 10 shows the reconstructed FMT image sequence through all time points ( 1 2 4 t t − , with a time interval 2.5 s) during FMT imaging process (0-60s). Furthermore, the second row of Fig. 10 depicts the time series of 3-D renderings of the reconstructed FMT image sequence. The entire tomographic image sequence (24 frame) is obtained by the proposed EKF-SOI reconstruction method. Fig. 10. The reconstructed FMT image sequence (total in 24 frames) from the phantom experiments based on one circle of projection data (24 projection data). The first shows the reconstructed 2-D tomographic sequence at different time points (2.5s, 15s, 30s, 45s, and 60s). The second row shows the 3-D rendering of the reconstructed images at the above time points. All images are displayed in the same range, respectively. Note that by the proposed EKF-SOI method, we can obtain entire tomographic images through all time points during FMT imaging processes.
The experimental results indicate that when using the conventional stationary method, the reconstructed image contain a high level of artifact [see Fig. 9(a)]. In contrast, when the proposed EKF-SOI reconstruction method is used, we can accurately resolve the ICG distribution in the phantom [see Fig. 9(b)]. In addition, it can also be found that compared with the EKF-SOI method, the reconstruction result from the EKF-DIM method contain more artifacts [see Fig. 9(c)]. On the other hand, similar to the numerical simulation, if using the conventional reconstruction method, only one tomographic image can be obtained based on the acquired 24 projection data. For comparison, based on the same data set, if using the EKF-SOI reconstruction method, we can acquire the entire 24 tomographic images corresponding to 1-24 angle during the imaging processes (see Fig. 10).  Figure 11 compares the reconstruction results from the in vivo experiment, obtained by the stationary, EKF-SOI, and EKF-DIM methods, respectively. Figure 11(a) shows the reconstruction image obtained by the stationary method, where the first circle of projection data was used as the input data of the reconstruction. Figure 11(b) shows the reconstruction results obtained by the proposed EKF-SOI method, which was imaged at 60s (corresponding to the 24th projection angle). For comparison, Fig. 11(c) shows the reconstruction results obtained by the EKF-DIM method at 60s.

In vivo experiment
For the in vivo reconstruction, the volume of interest was a 2.0 cm 2.2 cm 3.1cm × × 3-D region and sampled to 21 23 32 × × voxels. 5,132 voxels inside the imaged object and 5,117 source-detector pairs were used for reconstruction. In addition, to improve the reconstruction quality in the in vivo experiment, the reconstruction was performed by incorporating a heterogeneous forward model. Here, the forward model was constructed by assigning the midrange average optical properties to the heart, lungs, liver, and bones, which could provide an acceptable error of reconstructions in the mouse chest [26]. Similar to the simulation and phantom study, for the conventional stationary method, the reconstruction was performed by ART method. In the EKF reconstruction, we used 0 ρ 0 = , 0 0.025 = C I, 0 0.01 = Q I, and 0 0.0001 = R I, respectively. Furthermore, in the EKF-SOI reconstruction, the number of the second-order iteration was fixed to 5. These parameters were selected experimentally based on the results. The FMT reconstruction results obtained by the conventional methods, which are imaged at 0~1 min, 1~2 min, and 2~3 min after the injection of ICG, respectively. (e)-(i) The FMT reconstructed results obtained by the proposed EKF-SOI method based on the first projection data, which are imaged at 2.5s, 5.0s, 15s, 17.5s, and 30s, respectively, after the injection of ICG. (j)-(m) The EKF-SOI reconstructed results from the second circle, which are imaged at 62.5s, 65s, 92.5s, and 120s, respectively. Furthermore, Fig. 12 depicts the reconstruction images that are imaged after injecting ICG. Figures 12(b)-12(d) show the reconstruction results obtained by the conventional methods, which were imaged at 0~1 min, 1~2 min, and 2~3 min after the injection of ICG, respectively. Figures 12(e)-12(i) show the reconstructed result obtained by the proposed EKF-SOI method based on the first circle projection data, which were imaged at 2.5s, 5.0s, 15s, 17.5s, and 30s after the injection of ICG. Similarly, Figs. 12(j)-12(m) show the EKF-SOI reconstruction results from the second circle, imaged at 62.5s, 65s, 92.5s, and 120s, respectively. Fig. 13. The 3-D rendering of the reconstructed images at different time points (2.5s, 5s, 45s, 60s, and 62.5s) from the in vivo experiment. These reconstructed images are obtained by the proposed EKF-SOI method. All images are displayed in the same range.
Finally, Fig. 13 depicts the time series of 3-D FMT renderings acquired after the injection of ICG, which are obtained by the proposed EKF-SOI method. The reconstruction images at different time points (2.5s, 5s, 45s, 60s, and 62.5s) are selected to show the metabolic processes of ICG in the mouse liver.
By comparing with the anatomical structure of the liver provided by the XCT image, the experimental results indicate that the proposed EKF-SOI method can provide a good reconstruction quality [see Fig. 11(b)]. However, there is significant difference between the reconstruction results obtained by the EKF-SOI and EKF-DIM methods [see Fig. 11(b) and 11(c)]. Here, the reconstruction performance of EKF-DIM is significantly affected. On the other hand, similarly to the simulation and phantom experiments, if using the EKF-SOI method, based on one circle of projection data, we can obtain the entire 24-frame tomographic images corresponding to different projection angles (time points) (see Fig. 12 and Fig. 13). In contrast, if using the stationary method, based on the same data (one circle of projection data), only one tomographic image can be obtained [see Fig. 11(a)].

Discussion and conclusion
Dynamic FMT imaging allows fast and non-invasively resolving the 3-D distribution of fluorescence bio-markers within small animal in vivo, which is helpful for bio-medical research. However, reconstruction of the dynamic FMT is a challenging problem. The main reason is that the current reconstruction methods are performed based on the implicit assumption that the fluorophore's distributions are stationary during FMT imaging processes. To overcome the limitation, the main aim of this paper is to solve the non-stationary (timevarying) fluorescence tomography problem, which is achieved by extended kalman filter combined with second-order iterative method.
It could be observed from the numerical simulation, phantom, and in vivo experimental results that when the EKF-SOI method was used in reconstruction processes, we could clearly resolve the dynamic distributions of ICG in the imaged object (Figs. 5 and 9) compared with the conventional stationary method. In addition, we could also find that by the proposed EKF-SOI method, we could obtain the entire tomographic images at all time points based on one circle of projection data (Figs. 6, 10, and 13). For comparison, when using the conventional reconstruction method, only one tomographic image could be obtained based on the same projection data. That was, by using the proposed method, the reconstruction time resolution of FMT was reduced from 1 min (the stationary reconstruction method) to ~2.5 sec. As a result, by the proposed EKF-SOI method, it was possible to resolve the movement of fluorescence bio-markers within imaged object (Figs. 7 and 8) during the course of collecting one circle of projection data.
Based on the above results, we believe that the proposed EKF-SOI reconstruction method provides an attractive technique for solving the dynamic FMT problem. Firstly, this technique allows to achieve the entire tomographic images at all time points during FMT imaging process, which is helpful for imaging fast biological activities. In addition, the reconstruction time resolution of FMT only depends on the rotation time of the rotation stage rather than the data acquisition time for scanning one circle projection data. As a result, the imaging time resolution can be further reduced when a fast rotation speed is used. Furthermore, this method is fairly flexible to extend to other optical imaging modalities, e.g., photoacoustic tomography (PAT) [27] and x-ray luminescence computed tomography (XLCT) [28].
It should be noted that the reconstructed images at each time point obtained by EKF-SOI are based on the measurements corresponding to one angle projection data. It seems not obtain the reconstruction images from all angles, which may lead to a compromised spatial resolution and highly ill-posed inverse problem due to the limited projections. In addition, the initialization value of the EKF parameters (e.g., the fluorophore ρ , the estimation error covariance matrix t C , the noise covariance matrix t Q , and the covariance matrix of the observation errors t R ), are important, which may affect the performance of the proposed method. However, in this paper, the optimal EKF parameters are not determined. For example, in this work, the same value of 0 R is used for both the numerical simulation and real experiments. Considering that the measurement error of the real experiment is expected to be larger than the numerical simulation, a larger 0 R value might lead to better results for the EKF-SOI method in the real experiments. Moreover, the computation of the kalman smoother requires the all estimates from the forward run. Together, the forward matrix of FMT is generally large. As a result, the processing procedure of EKF-SOI is time-consuming, which is hard to meet the demands of dynamic imaging. Furthermore, to simplify the experiments, in the in vivo experiment, a thorough investigation of the reconstruction accuracy is not yet conducted. So, it is difficult to evaluate the performance of the proposed method in the in vivo experiment. To validate the reconstruction results from in vivo experiment, the multispectral epi-illumination cryoslicing imaging experiment may be performed [29]. These related work will be investigated in our future study.
In conclusion, we have applied the EKF-SOI method to FMT reconstruction to solve the time-varying FMT problem. The simulation, phantom, and in vivo experimental results indicate that based on the proposed method, we cannot only resolve the metabolic processes of fluorophore within the imaged object, but also detect the movement of fluorophore during the fast imaging processes, which will be helpful for resolving fast kinetic behaviors within small animals in vivo.