Open Access
Add to BibSonomy
Abstract
To alleviate the ill-posedness of diffuse fluorescence tomography (DFT) reconstruction and improve imaging quality and speed, a model-derived deep-learning method is proposed by combining extended Kalman filtering (EKF) with a long short term memory (LSTM) neural network, where the iterative process parameters acquired by implementing semi-iteration EKF (SEKF) served as inputs to the LSTM neural network correction model for predicting the optimal fluorescence distributions. To verify the effectiveness of the SEKF-LSTM algorithm, a series of numerical simulations, phantom and in vivo experiments are conducted, and the experimental results are quantitatively evaluated and compared with the traditional EKF algorithm. The simulation experimental results show that the proposed new algorithm can effectively improve the reconstructed image quality and reconstruction speed. Importantly, the LSTM correction model trained by the simulation data also obtains satisfactory results in the experimental data, suggesting that the SEKF-LSTM algorithm possesses strong generalization ability and great potential for practical applications.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Diffuse fluorescence tomography (DFT) with the merits of high sensitivity, non-invasiveness and non-radiation [1] demonstrates great potential for tumor diagnosis [2] drug development [3] and treatment evaluation [4]. However, due to the strong scattering of photons in biological tissue and incomplete surface measurements, DFT reconstruction is an ill-posed inverse problem [5], and the reconstructed images suffer from low spatial resolution, accuracy, as well as susceptibility to noise and model errors.
To alleviate the ill-posedness of DFT inverse problem, various strategies have been explored in an attempt to improve the quality of DFT reconstruction. The most commonly used strategy is known as regularization including Tikhonov regularization (L2 norm), Lp (0 < p ≤ 1) norm and total variation applied independently or jointly. Tikhonov regularization [6,7] is one of the most popular approaches to solve the inverse problem of DFT, in which a L2-norm constraint term is added into the data-fitting term to improve the stability, while it tends to cause artifacts and over-smoothed image boundary. As a sparsity constraint, L1 regularization [8,9] using prior knowledge of the sparse distribution of fluorescent sources is another effective strategy. However, there are some difficulties, such as over-sparseness, incomplete reconstruction of the fluorescent target, and lack of detailed information on the boundary. Apart from these, the quality of the DFT reconstruction can be improved by incorporating structural prior information obtained from XCT or MRI to the DFT reconstruction process, such as Helmholtz regularization [10], weighed segments regularization [11], and Laplace regularization [12].
Nevertheless, the above-mentioned methods usually lack a basis to adequately tackle the model mismatch due to not considering the measurement errors. In contrast, non-linear filtering scheme can handle measurement noise and also account for inaccuracies in forward models through an appropriate assignment of a process noise, and thereby provides a more rigorous framework to obtain estimation and error properties [13]. In diffuse optical tomography (DOT) field, Raveendran [14] et al applied ensemble Kalman filtering to retrieve DOT image by introducing the concept of pseudo-time, where the iterative solving process of image reconstruction was regarded as a dynamic changing process. The numerical simulation results showed that the position and shape of the targets could be reconstructed more accurately, compared with the regularized Gauss-Newton algorithm.
Baez [15] et al employed extended Kalman filtering (EKF) algorithm to estimate DOT image, and the results showed that the EKF-based method can effectively improve noise robustness compared with the NIRFAST solver (LM algorithm). In our group, Zhang [16] et al proposed a regularization-based EKF algorithm for DOT reconstruction, where Tikhonov regularization was incorporated to EKF process, resulting in superior imaging accuracy and noise robustness, especially under the circumstance of low absorption contrast and high noise level, compared with the conventional algebraic reconstruction technique (ART) and L2 regularization. Briefly, the strength of nonlinear filtering method is that it takes into account the influence of measurement noise on the filtering results during the iteration process, and thus the accuracy and noise robustness of image reconstruction can be improved to some extent. However, a primary disadvantage is that EKF tends to be time consuming due to multiple iterative updates and matrix inversions. In addition, similar to other traditional reconstruction strategies, EKF-based algorithm cannot substantially improve imaging quality, due to the inherent ill-posed nature of inverse problem.
Nowadays, neural networks have been introduced to the field of optical tomography reconstruction. Generally, the neural network-based methods can be categorized into end-to-end network and post-processing methods. The former directly establishes a nonlinear mapping relation from surface measurements to interior fluorescence distributions. For instance, Wang [17] et al proposed a stacked auto-encoder neural network to retrieve DFT image, and the numerical results demonstrated that the positions and shapes of the targets can be retrieved more accurately than the traditional ART algorithm. Li [18] et al proposed a framework of graph convolution network with ResNet-like architectures to achieve DFT image reconstruction with fewer parameters and higher speed. Zhang [19] et al proposed a 3D fusion dual-sampling deep learning network model in DFT to achieve ultra-high spatial resolution reconstruction. Though end-to-end network can enormously alleviate the ill-posedness of the inverse problem and improve image reconstruction quality, it turns the image reconstruction into a purely data-driven process, heavily relying on huge data to train the black box, and departing the physical relationship between measurement and reconstruction [20]. The alternative is post-processing method that generally first utilizes a conventional algorithm to obtain low-quality image, and then leverages neural network learning to correct the initial reconstruction image. Moreover, post-processing method can typically learn their mappings from less data, compared with purely data-driven model requiring a massive data set to learn a desirable mapping. Long [20] proposed a two-stage fluorescence tomographic reconstruction algorithm, where Tikhonov regularization was employed to obtain preliminary fluorophore distributions, sequentially the distributions of fluorophores were refined by convolutional neural networks (CNN). The results showed that the proposed method achieved more depth localization and less ambiguous reconstruction than L1 regularization. In photoacoustic tomography (PAT), Antholzer [21] et al applied filtered back projection algorithm to yield a low-quality image containing severe under-sampling artefacts, then CNN network was employed to map the intermediate reconstruction to an artefact-free final image. The results demonstrated that the proposed approach reconstructed images with a quality comparable to state-of-the-art iterative approaches for PAT from sparse data. Currently, the research on post-processing algorithms in the field of DFT is limited, and most studies are mainly in simulation stage.
As we mentioned above, EKF has been exploited for DOT image reconstruction and demonstrated its superiority over traditional deterministic methods. However, to our best knowledge, EKF has not been introduced to reconstruct DFT image. Herein, a model-derived deep-learning method, specifically, a post-processing method was proposed, where a semi-iteration EKF-based reconstruction was implemented to obtain the iterative process parameters, then LSTM neural network correct model was further employed to predict the optimal fluorescence distribution. It harnessed both the merits of EKF incorporating prior information and measurement errors to the model as well as long short term memory (LSTM) network with the good ability to mine crucial information from time series data via network learning. To verify the effectiveness of the proposed SEKF-LSTM algorithm for DFT imaging, a series of numerical simulations were conducted firstly, then the well-trained network model was applied to phantom and in vivo experiments, and the experimental results were quantitatively evaluated and compared with the EKF and semi-iteration EKF algorithms.
The rest of this work is organized as follows. In Section 2, the mathematical framework of DFT, EKF, and LSTM correction model are presented, and four metrics are introduced to quantitively assess the image quality. The methods and results of the simulation phantom and in vivo experiments are presented in Section 3 and Section 4, respectively. Finally, Section 5gives the conclusions.
2. Method
2.1 DFT technique
2.1.1 Forward problem
For continuous-wave DFT imaging, the light transport in biological tissue can be commonly described by using a set of the coupled diffusion equations as follows [22]:
Usually, the above equation can be resolved by combining with Robin boundary condition:
2.1.2 Inverse problem
To mitigate the influence of heterogeneous background optical properties as well as the errors between different measurement channels, normalized Born ratio method is used to reconstruct the fluorescence yields, which is written as follows:
Based on the FEM [23], Eq. (3) can be discretized into a matrix equation:
where $W({\boldsymbol r})$ is a M × N weight matrix with M and N representing the numbers of measurement data and finite element nodes, respectively; X represents the fluorescence yield $\eta {\mu _{af}}({\boldsymbol r})$ to be reconstructed.2.2 Extended Kalman filtering
2.2.1 State space model
The state space model of extended Kalman filtering consists of state and observation equations. In this work, the fluorescence yield is used as the state variable $X[k]$ to establish the state equation, and the ${I_{nb}}({\boldsymbol r})$ obtained from the measurement is applied to establish the observation equation, described by Eqs. (5) and (6), respectively.
Where $S\sim N(0,Q)$ and $D\sim N(0,R)$ denote the process noise and measurement noise, respectively, and both are zero-mean Gaussian white noise; Q and R represent the corresponding noise covariance, respectively. Herein, $W$ is referred to as measurement matrix that is consistent with the weight matrix in Eq. (4).2.2.2 DFT reconstruction based on an EKF model
The EKF solution to the above equation is given by a prediction-update procedure. In the prediction stage, the fluorescence yield ${\hat{X}^ - }[k|k - 1]$ and error covariance matrix ${\hat{P}^ - }[k|k - 1]$ at the kth step are predicted based on the one step ahead predictions $\hat{X}[k - 1]$ and $\hat{P}[k - 1]$, respectively. The estimation process can be formulated as follows:
where E is an identity matrix.In the update stage, the Kalman gain $G[k]$ at the kth step is calculated based on ${\hat{P}^ - }[k|k - 1]$ and measurement matrix W, as shown in Eq. (9). The fluorescence yield $\hat{X}[k]$ for step k is updated using the predicted ${\hat{X}^ - }[k|k - 1]$, $G[k]$, and the difference $\alpha [k]$ between the observed value and the predicted value, as shown in Eq. (11). The error covariance matrix $\hat{P}[k]$ for step k is updated using the predicted ${\hat{P}^ - }[k|k - 1]$, $G[k]$, and W, as shown in Eq. (12).
The above equations demonstrate that the fluorescence yield predicted at the next step is associated with the updated fluorescence yield ${\hat{X}^ - }[k|k - 1]$, the Kalman gain $G[k]$ and $\alpha [k]$ at the current step. Therefore, the three factors are used as the inputs to the LSTM correction model, and the output is the actual fluorescence yield $Y[k]$, as illustrated in Fig. 1. Accordingly, a mapping between the input and output can be expressed as:
Usually, the EKF parameters are empirically determined. In this work, we select $\boldsymbol X[0] = {10^{ - 6}}{\boldsymbol I_{N \times 1}}$, $P[0] = {10^{ - 8}}{\boldsymbol I_{N \times N}}$, $\boldsymbol Q = {10^{ - 8}}{\boldsymbol I_{N \times N}}$, and $\boldsymbol R = {10^{ - 9}}{\boldsymbol I_{1 \times 1}}$.
2.3 LSTM correction model
The process parameters iteratively obtained by the SEKF algorithm can be considered as time evolving parameters. Therefore, we selected the network that is suitable for time-series prediction. For learning temporal features in deep learning models, recurrent neural networks (RNN) have been a crucial approach for solving time series tasks. The long short-term memory (LSTM) network is a special form of recurrent neural network (RNN) that is capable of learning long term dependencies in data. The typical feature of RNN architecture is a cyclic connection, which enables the RNN to possess the capacity to update the current state based on past state and current input data, as shown in Fig. 2.
The hidden state ${h_{{\kern 1pt} t}}$ can be expressed as follows:
where ${W_x}$ and ${W_h}$ are the weights of the input ${x_t}$ and the hidden state of the last time ${h_{t - 1}}$, respectively, ${b_h}$ is bias, and $\tanh$ denotes hyperbolic tangent activation function.However, RNN is prone to suffer from gradient vanishing or exploding problems. Therefore, LSTM neural network is designed and applied [24].The LSTM structure leverages three kinds of multi-threshold gate units including forget gate ${f_t}$, input gate ${i_t}$, and output gate ${o_t}$, to suppress vanishing gradient and exploding gradient, as illustrated in Fig. 3. ${f_t}$ controls the amount of information that needs to be discarded when the network updates the state next time, ${i_t}$ controls the information flow of candidate information ${\tilde{c}_t}$ to the current cell state ${c_t}$ using a point-wise multiplication operation of ‘sigmoid’ and ‘tanh’, and ${o_t}$ decides which information should be passed on to the next hidden state.
LSTM calculation can be described by Eqs. (15)–(20):
where σ is sigmoid activation function.In this work, the LSTM network contains two LSTM layers and two fully connected (FC) layers, as shown in Fig. 4. The number of input neurons is determined by the EKF process parameters ${\hat{X}^ - }[k|k - 1]$, $G[k]$, and $\alpha [k]$; the number of output neurons is equal to that of finite element nodes; the number of hidden layer neurons in the LSTM and fully connected layers are set to 369 and 269, respectively. In addition, the dropout layer with a probability of 0.2 is added after the second LSTM layer and the first FC layer to prevent the model from over-fitting.
Mean square error (MSE) as the loss function is utilized to train network as follows:
2.4 Evaluation metrics
In order to quantitatively evaluate the performance of the EKF, SEKF and SEKF-LSTM algorithms, four evaluation metrics, namely, contrast-to-noise ratio (CNR), full width at half maximum (FWHM), quantitative reconstruction rate (QR), and mean absolute error (MAE) are adopted.
Roughly, QR represents quantitative accuracy of the reconstructed targets, FWHM implies the reconstruction spatial resolution, CNR is defined as a criterion for evaluating the contrast between the reconstructed target and background, i.e., the target visuality, and MAE demonstrates the accuracy of the whole reconstructed image. In general, better image reconstruction can be identified by a QR value closer to 1, a FWHM value closer to the true value, a larger CNR and a smaller MAE.
3. Experimental methods
To evaluate the performance of the proposed SEKF-LSTM algorithm, a series of numerical simulations, phantom and in vivo experiments were conducted and quantitatively compared with the semi-iteration EKF (SEKF) and full-iteration EKF- (abbreviated to EKF) algorithms.
3.1 Simulation experiments
3.1.1 Training dataset settings
Figure 5 shows a 2-D circular numerical phantom model with a radius of $R = 15mm$, embedded a circular target with a radius r. 16 source-detection pairs are evenly distributed around the phantom. The 16 sources are divided into four groups. When 4 modulated sources (NO.1-4) illuminate the phantom surface simultaneously, the remaining 12 fibers act as detection fibers to collect the transmission light. Sequentially, the phantom is excited by NO.5-8, NO.9-12 and NO.13-16 light sources, respectively, and the other fibers are taken as detection fibers accordingly. Consequently, a complete measurement can acquire a total of 16 × 12 data. The phantom is discretized into 721 finite element nodes and 1350 triangular elements. The reduced scattering coefficient and absorption coefficient of both the background and the target are set to ${u_s}^{\prime} = 1m{m^{ - 1}}$ ${u_a} = 0.0035m{m^{ - 1}}$, respectively. It is worth noting that the same experimental setup, including the geometry and optical coefficients, are also used in the subsequent phantom experiments, facilitating the mutual collation between the numerical and phantom experiments.
To generate simulation dataset, the fluorescence yields of the targets were varied from $0.0005\textrm{ }m{m^{ - 1}} - 0.003m{m^{ - 1}}$ in steps of $0.0005\textrm{ }m{m^{ - 1}}$, and the target with the radius of 2.5, 3.5, 4.5, or $5.5mm$ was randomly assigned 150 different positions over the phantom, resulting in a total of 3600 samples. To more closely simulate the real measurement, 40 dB Gaussian noise was added to the pure simulated data.
EKF iteration procedure includes inner-outer loops, where the number of inner iterations is equal to that of the measurement data, while the number of outer iterations is usually determined empirically. To further generate SEKF-LSTM network training samples, the semi-iteration EKF (SEKF) that stops iterating after executing two times of outer iterations was firstly implemented to obtain the process parameters (${\hat{X}^ - }[k|k - 1]$, $G[k]$, and $\alpha [k]$) from the 3600 sets of simulated data. Subsequently, in order to reduce model training time while maintaining model performance, the process parameters of the last ten iterations were selected as inputs to train the network, so the total number of samples for the LSTM correction model was 3600 × 10. The samples were split into training and validation dataset with a proportion of 7:3. The model performed 600 epochs of training using Adam algorithm with a batch size of 64 to complete the gradient optimization, and the learning rate was set to a fixed value of 0.0005. All computer processing was accomplished by a personal computer with an Intel Core i5-8250U CPU @ 1.60 GHz 1.80 GHz processor and 8GB of memory, and the time of training to convergence was about 1.7 h (Windows 10 and pytorch3.8).
3.1.2 Simulation scenarios settings
To comprehensively assess these algorithms, two main scenarios were considered:
- (1) Different target sizes: in order to probe the performance of the algorithms for reconstructing different target sizes, the radii of the target were determined to be 2, 3, and $4mm$, respectively; the fluorescence yield value of the target was fixed to $0.003m{m^{ - 1}}$.
- (2) Different target fluorescence yield values: in order to probe the performance of the algorithms for reconstructing different target fluorescence yields, the fluorescence yield values of the target were determined to be 0.001, 0.002 and $0.003m{m^{ - 1}}$, respectively; the radius of the target was fixed to $3mm$.
3.2 Phantom experiments
In order to further evaluate the performance of the three algorithms, phantom experiments were carried out based on our home-made DFT imaging system. In this system, lock-in photon-counting technique were employed to obtain high sensitivity, large dynamic range, and high ability to reject ambient light. The light sources are laser diodes (L785P25, Thorlabs) with a wavelength of 785 nm; the detectors are photomultiplier tubes (H8259-02e, Hamamatsu, Japan) having a limited dynamic range of 100 dB, and a complete measurement time for one frame image was 10.08 s. The technical details of the system can be referred to our previous work [25]. Figure 6(a) and (b) illustrate the sketch and the photo of the phantom, respectively. Herein, a custom-made cylindrical phantom made of polyoxymethylene with a diameter of $30mm$ and a height of $50mm$ was employed, where three circular holes of $4mm$ diameter and $35mm$ depth were drilled at the positions that were $9mm$ away from the central axis of the phantom and 120° separating angles among them. To mimic a tumor in biological tissue, the solution of 1% Intralipid mixed with indocyanine green (ICG, Liaoning Tianyi Biological Pharmaceutical Co, Ltd, China.) of 1.0, 3.0, and $5.0\mu mol/L$, respectively, was filled in No. 2 hole, while the other two holes (No.1 and 3) were filled with 1% Intralipid to match with the background. The quantum efficiency of ICG is 0.016, the peak excitation and emission wavelengths are 785 nm and 820 nm, respectively [26].
3.3 In vivo experiments
In order to verify the feasibility of the proposed algorithms in practical application, an in vivo experiment was conducted on an 8-week-old kunming mouse (weighting about 25 g) purchased from Tianjin Aoxin Experimental Animal Marketing Co. LTD (Tianjin, China). The animal experiment has been reviewed by the Animal Ethics Committee, and the procedures were developed in accordance with the guidelines of the Council for Control and Supervision of Animal Testing, Ministry of Health, PRC.
To perform the in vivo experiment, a transparent glass tube of $4mm$ inner diameter and $30mm$ height filled with ICG diluted with saline ($50\mu g/ml$) was employed as a fluorescent target. Before implementing the in vivo measurement, the mouse was anesthetized by intraperitoneal injection of 4% chloral hydrate at a dose of $0.01ml/g$, and the anesthetized mouse was depilated to reduce the effect of the hair for better detection of the excitation and emission light. Subsequently, the transparent glass tube was surgically buried into the mid-abdominal cavity of the mouse. The abdominal incision was closed with surgical sutures, and was applied to the suture to prevent incision infection. Finally, the mouse was placed in a polyformaldehyde cylinder chamber with a diameter of $23mm$, a height of $75mm$ and a wall thickness of $0.5mm$ surrounded by matching fluid (1% intralipid), as shown in Fig. 7(b). It is worth mentioning that the training set was regenerated by using the absorption coefficients of $0.035m{m^{ - 1}}$[27], corresponding to the average value of bulk soft tissues and the background radius of $R = 11.5mm$ to match the in vivo experiment, since the absorption coefficient of the in vivo mouse and the radius of the imaging chamber were different from those of the numerical phantom.
Fig. 7. (a) In vivo mouse surgically buried transparent glass tube in mid-abdominal cavity; (b) in vivo measurement with DFT imaging system.
4. Results
4.1 Simulation results
4.1.1 Different target sizes
Figure 8 shows the reconstructed images with the fixed target fluorescence yield value of $0.003m{m^{ - 1}}$, while variable target radii of 2, 3, and $4mm$, respectively. The profile plotted along the dashed line passing through the centers of the background and the target quantitatively depicts the fluorescence yield values. In order to facilitate the comparison of different algorithms, the images reconstructed under the same fluorescence yield value are shown using the same colormap. Thus, all the reconstructions are presented in the same color scale to allow level comparison in this case. Visually, we can observe that the reconstructed targets for all the algorithms are clearly distinguishable from the background and demonstrate strong noise robustness attributed to incorporating measurement noise and prior information in EKF-based reconstruction. With the increase of target size, the target becomes more distinguishable, suggesting larger fluorescence yield value. As expected, the profiles at the rightmost column demonstrate the SEKF achieves the lowest fluorescence yield due to insufficient convergence, especially in the case of small target size ($r = 2mm$). In contrast, the fluorescence yield values obtained by the SEKF-LSTM are closer to the true value under different target sizes. From the view of the target shape and size, the SEKF-LSTM achieves sharper edge and more faithfully shape, compared with the SEKF and EKF algorithms suffering from over-smooth and a certain deformation of the target shape.
Fig. 8. The reconstructed images and profiles of fluorescent yields with the variable target radius of 2, 3, and 4 mm, and the fixed fluorescence yield of 0.003mm-1, obtained by EKF, SEKF, and SEKF-LSTM, respectively.
Furtherly, the performances of different algorithms were quantitatively evaluated by the metrics of CNR, FWHM, QR and MAE according to Eqs. 22–24, and the column diagrams of the metrics were plotted. As shown in Fig. 9, the SEKF-LSTM algorithm achieves distinctly high CNRs for different target sizes, indicating the SEKF-LSTM gains the optimal noise robustness. With the increase of the target size, the FWHMs of the profiles obtained by all the algorithms are closer to the true ones. Nevertheless, the traditional EKF-based algorithms, especially the SEKF, demonstrate an overestimation of the target size, in comparison to the SEKF-LSTM algorithm where the FMHMs are approximately equal to the true values for all the cases. As expected, the QR of the SEKF are lower than that of the EKF due to incomplete convergence and demonstrate an increasing tendency with the increasing of target size. Comparatively, the QR of the SEKF-LSTM algorithm demonstrates a more stable performance (close to 1) under different target size. In addition, the MAEs obtained by the SEKF-LSTM are far less than those of the EKF-based algorithm, indicating that the reconstructed images are more faithful the true ones. Comprehensively, the SEKF-LSTM algorithm possesses higher noise robustness, spatial resolution, quantitativeness and accuracy.
4.1.2 Different target fluorescence yields
Figure 10 illustrates the reconstructed images and the profiles for targets with the fixed radius of $3mm$, while variable fluorescence yield values of 0.001, 0.002, and $0.003m{m^{ - 1}}$, respectively. Visually, all the algorithms can successfully estimate the target location and shape, and the SEKF-LSTM-based algorithm achieve the sharpest edge and the most faithful size in target reconstruction. From the profiles at the bottom of the Fig. 10, we can see that the reconstructed fluorescence yields demonstrate an increasing trend with the increase of the fluorescence yield values for all the algorithms. Furthermore, the profiles at the rightmost column display that the fluorescence yield values obtained by SEKF-LSTM algorithm are the closest to the true ones and present the highest quantitativeness under different fluorescence yields. However, a certain degree of overestimation exists in the scenario of low fluorescence yield value of $0.001m{m^{ - 1}}$, in comparison to the other two algorithms demonstrating obvious underestimation.
Fig. 10. The reconstructed images and profiles of fluorescent yield with the variable target fluorescence yields of 0.001, 0.002, and 0.003mm-1, and the fixed radius of 3 mm, obtained by EKF, SEKF, and SEKF-LSTM respectively.
Correspondingly, Fig. 11 displays the column diagrams of the four metrics at different fluorescence yields of 0.001, 0.002, and $0.003m{m^{ - 1}}$, obtained using the EKF, SEKF and SEKF-LSTM algorithms, respectively. It can be seen that the SEKF-LSTM algorithm obtains the highest CNRs, and the FWHMs are the closest to the trues ones. In comparison to the QRs below 0.8 and 0.6 obtained by the EKF and SEKF algorithms respectively, all the QRs obtained by the SEKF-LSTM algorithm are larger than 0.8. Furtherly, we can see that the MAEs obtained by all algorithms show a rise trend with the increasing of fluorescence yield values. Nevertheless, the MAEs obtained by the SEKF-LSTM algorithm are generally one order of magnitude smaller than the EKF algorithms, indicating that the proposed method has a superior performance on fidelity. These assessment results jointly verify that our proposed SEKF-LSTM algorithm effectively improved the DFT image quality.
Fig. 11. The column diagrams of the CNRs, FWHMs, QRs and MAEs of the reconstructed images in Fig. 10.
Finally, the reconstruction times of the three algorithms were compared. To obtain the DFT image, the EKF-based algorithm took about 31.2144s, while the SEKF-LSTM algorithm took about 10.1033s + 0.0078s, including the reconstruction times of SEKF (10.1033s) and LSTM (0.0078s). Therefore, the reconstruction times of SEKF-LSTM and SEKF are comparable and obviously less than that of EKF.
4.2 Phantom experiment results
Figure 12 illustrates the reconstructed fluorescence yield images for the targets with the ICG concentrations of 1.0, 3.0, and $5.0\mu mol/L$, respectively, and the profiles plotted along the centers of the backgrounds and targets. Similar to numerical results, all the algorithms can successfully distinguish the targets from the background under different ICG concentrations, and the estimated target location demonstrates a good agreement with the true one. The targets shape obtained by the traditional EKF-based algorithm appears obvious expansion and slight deformation, while the targets estimated by the SEKF-LSTM-based algorithm are consistent to the true ones. As can be seen from the rightmost profiles, the SEKF-LSTM-based algorithm achieves the largest fluorescence yield values under different ICG concentrations. From the bottom profiles, we can see that the fluorescence yields reconstructed by the three algorithms increases with the increase of ICG concentration, while the increase from $3.0\mu mol/L$ to $5.0\mu mol/L$ is not obvious, probably due to fluorescence quenching effect.
Fig. 12. The reconstructed images and profiles of fluorescent yield with the variable target ICG concentrations of 1.0, 3.0, and 5.0µmol/L obtained by EKF, SEKF, and SEKF-LSTM respectively.
Because the true fluorescence yield values are unknow, Fig. 13 only displays the column diagrams of the CNRs and FWHMs obtained using the EKF, SEKF and SEKF-LSTM algorithms, respectively. Consistent with the numerical results, the CNRs of the SEKF-LSTM algorithm are the largest, and the FMHWs are the closest to the true value over the whole range of ICG concentrations. Additionally, the CNRs obtained by the EKF and SEKF are comparable, while the FWHMs obtained by the SEKF are larger than that of EKF. Comprehensively, these phantom experimental results are consistent with the simulation ones.
4.3 In vivo experiment results
Figure 14 illustrates the reconstructed in vivo images and profiles using the EKF, SEKF, and SEKF-LSTM algorithms, respectively. It can be observed that the targets acquired by the three algorithms are clearly visible, and the reconstructed target based on SEKF-LSTM algorithm has more regular shape and clearer edges. The profiles demonstrate the fluorescence yield values achieved by the SEKF-EKF are the largest, followed by the EKF and SEKF. Figure 15 shows the CNRs and FWHMs of the reconstructed images in Fig. 14. As can been seen in Fig. 15, the CNR obtained by the EKF, SEKF, and SEKF-LSTM algorithms are 2.78, 2.97, and 14.87, respectively, indicating that the SEKF-LSTM algorithm has the best anti-noise ability. The corresponding FWHMs are 5.3, 5.4, and $3.8mm$. Notably, the target size estimated by the SEKF-LSTM is closer to the glass tube diameter of $4mm$, demonstrating high spatial resolution. Therefore, the in vivo experimental results verify the consistencies between in vivo experiment and the phantom experiment.
Fig. 14. The reconstructed images and profiles obtained by EKF, SEKF and SEKF-LSTM in the in vivo experiment, respectively.
5. Discussion and conclusion
In this work, a post-processing algorithm combining semi-iteration EKF and LSTM neural network was proposed to improve the quality of DFT image reconstruction and the speed of imaging. To verify the performance of the proposed algorithm, the evaluation was performed with simulation experiments firstly, and quantitatively compared with the semi-iteration EKF and full-iteration EKF-based algorithms. The numerical simulation results showed the SEKF-LSTM achieved the optimal performance in terms of noise robustness, spatial resolution, accuracy, and cost less computation time under different target size and fluorescence yield. Subsequently, the phantom data probed by our homemade measurement system was further evaluated relayed on the well-trained network models using numerical simulation dataset. The target size, shape, position estimated by the SEKF-LSTM algorithm were in good agreement with the true one, and achieved the higher quantitativeness and spatial resolution in the scenario of different ICG concentrations, confirming the strength of the proposed algorithm. Finally, the in vivo experimental results obtained by the SEKF-LSTM demonstrated good consistency with the phantom experiment, suggesting the great potential for clinical applications. In summary, comprehensive experimental results validated that the SEKF-LSTM-based DFT algorithm produced superior imaging quality. More importantly, this work demonstrates that the network trained on the simulated dataset is also suitable for experimental data, suggesting the trained LSTM network has good generalization ability and great potential for practical applications. Though the image reconstruction time was significantly reduced in the SEKF-LSTM algorithm, the network training process of LSTM is time consuming. Because LSTM is inherently recursive, i.e., it needs to compute the current state based on the previous state, it cannot be trained in parallel. Gate Recurrent Unit (GRU) is a simplified version of LSTM, which reduces the number of parameters, and thereby increases computational efficiency compared with the LSTM. In the future work, we will further explore the feasibility of SEKF in conjunction with GRU or other networks suitable for time-series data prediction to improve DFT image reconstruction.
Funding
National Natural Science Foundation of China (62075156, 62175183, 62205239, 81671728, 81971656, 82071971); Tianjin Municipal Education Commission (No. 2021KJ264).
Disclosures
The authors declare no financial or commercial conflict of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. P. Zhang, C. B. Ma, F. Song, et al., “A review of advances in imaging methodology in fluorescence molecular tomography,” Phys. Med. Biol. 67(10), 10TR01 (2022). [CrossRef]
2. L. M. Habimana-Griffin, D. Z. Ye, J. Carpenter, et al., “Intracranial glioma xenograft model rapidly reestablishes blood–brain barrier integrity for longitudinal imaging of tumor progression using fluorescence molecular tomography and contrast agents,” J. Biomed. Opt. 25(02), 1 (2020). [CrossRef]
3. K. D. Zhou, Y. C. Ding, I. Vuletic, et al., “In vivo long-term investigation of tumor bearing mKate2 by an in-house fluorescence molecular imaging system,” Biomed Eng Online 17(1), 187 (2018). [CrossRef]
4. V. Ntziachristos, E. A. Schellenberger, J. Ripoll, et al., “Visualization of antitumor treatment by means of fluorescence molecular tomography with an annexin V-Cy5.5 conjugate,” Proc. Natl. Acad. Sci. U.S.A. 101(33), 12294–12299 (2004). [CrossRef]
5. H. G. Yi, S. H. Ma, R. G. Yang, et al., “Neighbor-based adaptive sparsity orthogonal least square for fluorescence molecular tomography,” J. Biomed. Opt. 28(06), 066005 (2023). [CrossRef]
6. H. Zhang, H. B. Guo, S. C. Li, et al., “L1–L2 Minimization Via A Proximal Operator For Fluorescence Molecular Tomography,” 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), Mexico, 3640–3645 (2021).
7. Y. H. Liu, S. X. Jiang, J. Liu, et al., “Reconstruction method for fluorescence molecular tomography based on L1-norm primal accelerated proximal gradient,” J. Biomed. Opt. 23(08), 1–11 (2018). [CrossRef]
8. J. C. Baritaux, K. Hassler, and M. Unser, “An efficient numerical method for general Lp regularization in fluorescence molecular tomography,” IEEE Trans. Med. Imaging 29(4), 1075–1087 (2010). [CrossRef]
9. H. F. Wang, C. Bian, L. X. Kong, et al., “A Novel Adaptive Parameter Search Elastic Net Method for Fluorescent Molecular Tomography,” IEEE Trans. Med. Imaging 40(5), 1484–1498 (2021). [CrossRef]
10. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, et al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express 15(13), 8043–8058 (2007). [CrossRef]
11. D. Hyde, E. L. Miller, D. H. Brooks, et al., “Data specific spatially varying regularization for multimodal fluorescence molecular tomography,” IEEE Trans. Med. Imaging 29(2), 365–374 (2010). [CrossRef]
12. H. Meng, K. Wang, Y. Gao, et al., “Adaptive Gaussian Weighted Laplace Prior Regularization Enables Accurate Morphological Reconstruction in Fluorescence Molecular Tomography,” IEEE Trans. Med. Imaging 38(12), 2726–2734 (2019). [CrossRef]
13. M. Khodarahmi and V. Maihami, “A Review on Kalman Filter Models,” Archives of Computational Methods in Engineering 30(1), 727–747 (2023). [CrossRef]
14. T. Raveendran, S. Gupta, R. M. Vasu, et al., “A pseudo-time EnKF incorporating shape based reconstruction for diffuse optical tomography,” Med. Phys. 39(2), 1092–1101 (2012). [CrossRef]
15. G. R. Baez, J. A. Pomarico, and G. E. Elicabe, “An improved extended Kalman filter for diffuse optical tomography,” Biomed. Phys. Eng. Express 3(1), 015013 (2017). [CrossRef]
16. J. Zhang, L. M. Zhang, Z. Y. Liu, et al., “Tikhonov regularization-based extended Kalman filter technique for robust and accurate reconstruction in diffuse optical tomography,” J. Opt. Soc. Am. A 40(1), 10–20 (2023). [CrossRef]
17. H. Q. Wang, J. Y. Gao, Z. Zhao, et al., “Fluorescence Molecular Tomography Reconstruction of Small Targets Using Stacked Auto-Encoder Neural Networks,” IEEE Access 8, 37657–37663 (2020). [CrossRef]
18. D. S. Li, C. X. Chen, J. F. Li, et al., “Reconstruction of fluorescence molecular tomography based on graph convolution networks,” J. Opt. 22(4), 045602 (2020). [CrossRef]
19. P. Zhang, G. D. Fan, T. T. Xing, et al., “UHR-DeepFMT: Ultra-High Spatial Resolution Reconstruction of Fluorescence Molecular Tomography Based on 3-D Fusion Dual-Sampling Deep Neural Network,” IEEE Trans. Med. Imaging 40(11), 3217–3228 (2021). [CrossRef]
20. F. Long, “Deep learning-based mesoscopic fluorescence molecular tomography: an in silico study,” J. Med. Imag. 5(03), 1 (2018). [CrossRef]
21. S. Antholzer, M. Haltmeier, and J. Schwab, “Deep learning for photoacoustic tomography from sparse data,” Inverse Probl. Sci. Eng. 27(7), 987–1005 (2019). [CrossRef]
22. Y. Q. Zhang, L. M. Zhang, G. Y. Yin, et al., “Assessing indocyanine green pharmacokinetics in mouse liver with a dynamic diffuse fluorescence tomography system,” J. Biophotonics 11(10), e201800041 (2018). [CrossRef]
23. F. Gao, H. J. Zhao, Y. Tanikawa, et al., “A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography,” Opt. Express 14(16), 7109–7124 (2006). [CrossRef]
24. K. Smagulova and A. P. James, “A survey on LSTM memristive neural network architectures and applications,” Eur. Phys. J. Spec. Top. 228(10), 2313–2324 (2019). [CrossRef]
25. L. M. Zhang, N. Cheng, H. Liu, et al., “High-sensitivity dynamic diffuse fluorescence tomography system for fluorescence pharmacokinetics,” J. Biomed. Opt. 27(04), 046002 (2022). [CrossRef]
26. T. Desmettre, J. M. Devoisselle, and S. Mordon, “Fluorescence Properties and Metabolic Features of Indocyanine Green (ICG) as Related to Angiography,” Surv Ophthalmol 45(1), 15–27 (2000). [CrossRef]
27. D. Hyde, R. Schulz, D. Brooks, et al., “Performance dependence of hybrid x-ray computed tomography/fluorescence molecular tomography on the optical forward problem,” J. Opt. Soc. Am. A 26(4), 919–923 (2009). [CrossRef]
