Gradient Projection for Sparse Reconstruction Method for Dynamic Fluorescence Molecular Tomography

Abstract

Dynamic fluorescence molecular tomography (FMT) is a promising optical imaging technique for three-dimensionally demonstrating the metabolic process of fluorophore in small animals. Conventional FMT methods focus on reconstructing static distribution of fluorescent yield, and the reconstruction results may perform poorly if the boundary measurement data, acquired from time-varying fluorophore, were directly used in these methods. In this study, we apply joint $\ell_{1}$ and Laplacian manifold regularization model to dynamic FMT. The $\ell_{1}$-norm regularization method is used to deal with the ill-posedness of FMT, and the Laplacian manifold regularization is introduced to obtain spatial structure information of boundary measurements. Then, we use gradient projection for sparse reconstruction (GPSR) method to solve the joint regularization model. Since the boundary measurements are obtained from different time points, the input data is converted from a vector to a matrix, and each column of the matrix corresponds to a time point. Thus, a sequence of fluorophore concentration images, corresponding to different time points, can be reconstructed in one step. Numerical simulation experiments are performed and the results indicate that the proposed method can recover the dynamic fluorophore well.

1 Introduction

Dynamic fluorescence molecular tomography (FMT), as a promising optical tomographic imaging technique, provides richer information than conventional static imaging, and has the ability to capture the absorption, distribution, and elimination of fluorescent bio-markers within small animals [1,2,3]. It is helpful in better studying drug delivery and disease progression.

In static FMT, the traditional methods [4,5,6] assumed the concentration of fluorescent probe is constant, so these methods only reconstruct static distribution of fluorescent yield from data sets acquired within several minutes, but fail to capture the temporal information of the time-varying fluorophore. Besides, the reconstructed results perform poorly if the measurements, corresponding to the non-stationary fluorophore concentration, are directly used in the static FMT reconstruction, especially when the fluorophore concentrations vary rapidly.

To address the problem, some theoretical and experimental solutions have been proposed [7,8,9,10], and the fluorophore is assumed as a discrete time stochastic process. These methods allows the fluorophore concentration change over time, and they are capable to solve the dynamic FMT reconstruction problem. Since the fluorephore concentrations, viewed as a random variable at each time point, are temporally correlated, the acquired measurement data are highly correlated with one another along the time axis. These methods, however, obtain measurement data from different time first, and then reconstruct a series of static distributions of fluorescent yield respectively, fail to make full use of the temporal correlations of boundary measurements.

To overcome these limitations, in this paper, we propose a joint regularization model to reconstruct a sequence of FMT images of different time points in one step. Both accurate transportation model and sophisticated inverse algorithm are indispensable in FMT. To deal with the high ill-posedness of FMT, many regularization methods have been introduced to solve the inverse problem. $\ell_{2}$-norm regularization, such as Tikhonov regularization [11, 12], was widely used in FMT. However, $\ell_{2}$-norm regularization tends to incur over-smoothness in reconstructed image. To improve the imaging quality, $\ell_{1}$-norm regularization, as an alternative solution, was used to promote the sparsity of the solution [13,14,15]. For dynamic FMT, spatial structure is significant information for reconstructing the distribution of fluorophore. To utilize both sparsity and spatial structure information, manifold-based learning method [16,17,18,19] is needed. Thus, we introduce Laplacian manifold regularization and combine it with $\ell_{1}$-norm regularization to improve dynamic FMT reconstruction. Then, we use gradient projection for sparse reconstruction (GPSR) method [20, 21] to solve the joint regularization model.

The outline of this paper is as follows. In Sect. 2, the methods are detailed. Section 3 presents the experimental materials, and the reconstruction results are shown. In Sect. 4, we draw conclusions based on the simulation results.

2 Methods

2.1 Dynamic FMT Forward Model

For dynamic FMT, the photon propagation in biological tissues can be modeled using the diffusion equation (DE). In a continuous-wave domain with point excitation sources, the propagation of excitation and emission lights can be presented as follows:

$$ \left\{ {\begin{array}{*{20}c} {\nabla \cdot [D_{x} (r)\nabla\Phi (r,t)] - \mu_{ax} (r)\Phi _{x} (r,t) = -\Theta \delta (r - r_{s} ,t)} \\ {\nabla \cdot [D_{m} (r)\nabla\Phi _{m} (r,t)] - \mu_{am} (r)\Phi _{m} (r,t) = -\Phi _{x} \eta \mu_{\alpha f} (r,t)} \\ \end{array} } \right.\begin{array}{*{20}c} {(r \in\Omega )} \\ {(r \in\Omega )} \\ \end{array} $$

(1)

where $\Omega$ is the domain of the imaged object, $r$ is the position vector. $\mu_{\alpha x,m} (r)$ represents the absorption coefficient, $t$ denotes time. The subscript $x$ and $m$ represent the excitation process and emission progress respectively. $\Phi _{x,m} (r,t)$ denotes photon intensity and $D_{x,m} (r)$ is the diffusion coefficient of the tissue at position $r$. $\delta (r - r_{s} ,t)$ is the point excitation source and $\Theta$ is the amplitude. $r_{s}$ is the position vector of different excitation sources. $\eta$ denotes the fluorescent yield efficiency, $\mu_{\alpha f} (r,t)$ is the absorption coefficient of fluorescent to the excitation light.

In this paper, finite element method (FEM) is chosen to solve the diffusion equations that the imaged domain is discretized into a mesh with $N$ nodes. Equation (1) is converted into the following linear equations:

$$ \left\{ {\begin{array}{*{20}c} {K_{x}\Phi _{x} = S_{x} } \\ {K_{m}\Phi _{m} = FX} \\ \end{array} } \right. $$

(2)

where $K_{x}$ and $K_{m}$ represent the system matrix in excitation process and emission process, respectively. $S_{x}$ denotes the excitation source distribution, $F$ is an $N \times N$ symmetric matrix. Based on Eq. (2), the final equation can be formulated as:

$$ AX =\Phi $$

(3)

Equation (3) shows the linear relationship between the unknown fluorescent yield $X$ and the boundary measurements $\Phi$. $A$ is an $M \times N$ weight matrix depending on the geometry and optical parameters.

For static FMT, the fluorophore is assumed to be stationary that it does not change during FMT imaging processes. Thus, $X$ is a $N \times 1$ vector represents the fluorescent yield at a given time $t$, and $\Upphi$, also a $N \times 1$ vector, is the corresponding boundary measurements. However, in dynamic FMT, a time series of boundary measurements are acquired at intervals of minutes over time and the fluorescent yield can be denoted as $X = \{ X_{1} ,X_{2} , \cdots ,X_{t} , \cdots ,X_{k} \}$, where $k$ is the number of time points. As a consequence, both of $X$ and $\Phi$ are $N \times k$ matrices.

2.2 Dynamic FMT Inverse Model

The static reconstruction algorithms fail to solve the time-varying measurements, in this paper, we proposed a joint $\ell_{1}$ and Laplacian manifold regularization model to improve the reconstruction performance.

To cope with the high ill-posed nature of Eq. (3) and the existence of noise, some form of regularization method is indispensable in the inversion. Compared with $\ell_{2}$-norm regularization method, $\ell_{1}$-norm regularization is used in this paper due to its sparsity. The inverse problem can be formulated as:

$$ \arg \mathop {\hbox{min} }\limits_{X} \frac{1}{2}||AX -\Phi _{m} ||_{2}^{2} + \tau ||X||_{1} $$

(4)

In order to improve the quality of reconstruction, both structural priors and temporal correlations of boundary measurements should be fully utilized. Thus, we introduce the Laplacian manifold regularization into the inverse process:

$$ \arg \mathop {\hbox{min} }\limits_{X} \frac{1}{2}||AX -\Phi _{m} ||_{2}^{2} + \tau ||X||_{1} + \frac{\lambda }{2}X^{T} LX $$

(5)

where $X^{T} LX$ is the Laplacian manifold regularization term and $L$ is the regularization.

Let $p_{i} (i = 1,2, \ldots ,N)$ denotes the space nodes, $e_{i,j} (i,j = 1,2, \ldots ,N)$ denotes the edges of FEM mesh. $e_{i,j} = 1$ represents there is an edge between $p_{i}$ and $p_{j}$, otherwise, $e_{i,j} = 0$.

The Laplacian manifold regularizer is defined as:

$$ \sum\limits_{i = 1}^{N} {\sum\limits_{j = 1}^{N} {w_{i,j} } } (x_{i} - x_{j} )^{2} $$

(6)

where $w_{i,j}$ is the weight of the edge between node $p_{i}$ and $p_{j}$, defined as:

$$ w_{i,j} = \left\{ {\begin{array}{*{20}c} {\exp ( - \frac{{||p_{i} - p_{j} ||_{2}^{2} }}{{\sigma^{2} }})} & {{\text{if e}}_{i,j} = 1} \\ 0 & {{\text{if i = j or e}}_{i,j} = 0} \\ \end{array} } \right. $$

(7)

where $\sigma > 0$ is the parameter to adjust the weight matrix. It is apparent that $w_{i,j} \in [0,1]$, and the closer $p_{i}$ from $p_{j}$, the closer $w_{i,j}$ to 1. By simple mathematical derivation, we can obtain the formula:

$$ \sum\limits_{i = 1}^{N} {\sum\limits_{j = 1}^{N} {w_{i,j} } } (x_{i} - x_{j} )^{2} = 2X^{T} (D - W)X $$

(8)

where $W$ is a symmetrical weight matrix, which equals to $(w_{ij} )_{N \times N}$. $D$ is a diagonal matrix, which equals to $diag(d_{1} ,d_{2} , \ldots ,d_{N} )$. Let $L = D - W$, then we obtain the Laplacian manifold regularization item of Eq. (5).

2.3 Gradient Projection for Sparse Reconstruction Method

To resolve the dynamic optical imaging problem, the conventional approach estimates each image independently, and then assembles these results into a time sequence. The reconstruction results show good quality at each time point, but it fails to make full use of the measurements $\Phi = \{\Phi _{1} ,\Phi _{2} , \ldots ,\Phi _{k} \}$ that every column is temporarily correlated with one another. In this paper, GPSR is utilized to solve Eq. (5).

Through GPSR approach, the variable $X$ is split into positive and negative parts. Equation (5) can be rewritten as the following quadratic program:

$$ \arg \hbox{min} \frac{1}{2}Z^{T} BZ + C^{T} Z\,\,s.t.Z \ge 0 $$

(9)

where $Z = \left[ {\begin{array}{*{20}c} {Z_{p} } \\ {Z_{n} } \\ \end{array} } \right]$, $C = \left[ {\begin{array}{*{20}c} {\tau 1_{N} - A^{T}\Phi } \\ {\tau 1_{N} + A^{T}\Phi } \\ \end{array} } \right]$, and $B = \left[ {\begin{array}{*{20}c} {A^{T} A + \lambda L} & { - (A^{T} A + \lambda L)} \\ { - (A^{T} A + \lambda L)} & {A^{T} A + \lambda L} \\ \end{array} } \right]$.

These relationships are satisfied by $Z_{p} = (x_{i} )_{ + }$ and $Z_{n} = ( - x_{i} )_{ + }$ for all $i = 1,2, \ldots ,N$, where $( \cdot )_{ + }$ denotes the positive-part operator. Then gradient projection algorithm is used to solve the problem in Eq. (9).

Let $F(Z) = 1/2Z^{T} BZ + C^{T} Z$, we search from each iterate $Z^{(v)}$ along the negative gradient $- \nabla F(Z^{(v)} )$, and iterate $Z^{(v)}$ evolves to iterate $Z^{(v + 1)}$ as follows:

$$ Z^{(v + 1)} = (Z^{(v)} - \alpha^{(v)} \nabla F(Z^{(v)} ))_{ + } $$

(10)

where $\alpha^{(v)} > 0$ is the scalar parameter. We define the vector $g^{(v)}$ as:

$$ g_{i}^{(v)} = \left\{ {\begin{array}{*{20}c} {[\nabla F(Z^{(k)} )]_{i} } & {if\begin{array}{*{20}c} {z_{i}^{(v)} > 0} & {or\,[\nabla F(z^{(v)} )]_{i} < 0} \\ \end{array} } \\ 0 & {otherwise} \\ \end{array} } \right. $$

(11)

The initial value of $\alpha^{(v)}$ can be computed as:

$$ \alpha_{0} = \frac{{(g^{(v)} )^{T} g^{(v)} }}{{(g^{(v)} )^{T} Bg^{(v)} }} $$

(12)

Compared with other methods, the proposed method has several features: First, the joint $\ell_{1}$ and Laplacian manifold regularization model makes full use of the temporal correlations of boundary measurements. Second, the input data $\Phi$ is a $N \times 4$ matrix, each column corresponds to a set of boundary measurements that obtained from different time points. The proposed algorithm could output the reconstructed results $X$, a $N \times 4$ matrix, in one step.

2.4 Evaluation Index

To quantify the reconstruction performance, position error (PE) was calculated to measure the distance variation between the reconstructed region and real region. The definition of PE is:

$$ PE = ||P_{r} - P_{0} ||_{2} $$

(13)

where $P_{r}$ and $P_{0}$ are the center coordinates of reconstructed and actual fluorescent sources respectively. We also calculated the volume of the reconstructed fluorophore, to compare with the volume of real target. The results is shown in Sect. 3.

3 Experiments

3.1 Experimental Setup

Numerical simulation is implemented to validate the performance of the joint regularization model and the proposed algorithms. A virtual mouse atlas was employed to provide the 3-D information, the mouse torso from the neck to the base of the kidneys was selected as the investigated region. As shown in Fig. 1(a), the simulation model includes five organs: muscle, heart, lungs, liver and kidneys. To simplify the simulation of the metabolic processes of drug, we focus on imaging kinetic behavior in liver.

Fig. 1.

(a) The mouse 3-D geometry model used in simulation studies. (b) Shows the position of the spherical target.

The spherical target is set in the liver and the center of it is at (16.4, 11, 16.4) mm, as shown in Fig. 1(b). The positions of four excitation point sources is located at Z = 16 mm. For each excitation location, the fluorescence is measured from the opposite side within 120 deg field of view.

The absorption and the scattering coefficients are essential in the photon propagation simulation. The following Table 1 gives the optical parameters of the mouse organs. The optical properties outside these organs were regarded as homogeneous.

Table 1. Optical parameters of the mouse organs.

To simulate the metabolic process, a series of spherical targets with different concentrations and sizes should be set as different fluorescents. We planned to set three simulation experiments: (a) Spherical targets with same concentration and different radii. (b) Spherical targets with different concentrations and same radius. (c) Spherical targets with different concentrations and different radii.

In this paper, we mainly focus on testing the proposed model and the algorithm, only the first experiment is illustrated, the other experiments would be done in the future work. According to the ICG concentration curves in [22], which mimic the metabolic processes of ICG in different organs and tissues, we set four spherical targets with four different radii: 1.0, 1.5, 2.0, 2.5 (mm) to simulate the diffusion progress of fluorescent at the four time points: 10, 15, 30, 60 (min). The variation of target concentration within each time point is neglected for simplification.

3.2 Forward Simulation

The photo propagation model of dynamic FMT reconstruction is illustrated in this section. Four spherical targets with different radii: 1.0, 1.5, 2.0, 2.5 (mm) and same concentration was set in the liver, as shown in Fig. 2.

Fig. 2.

Four mouse models with different spherical targets. The radius are 1.0, 1.5, 2.0 and 2.5 mm, corresponding to (a) to (d), to simulate the diffusion progress of fluorescent.

The photon density on the mouse surface was simulated using the FEM approach. In order to obtain the measurements, the four models was discretized into 10689, 10793, 10869, 10750 nodes, respectively. Figure 3 shows the forward results. 3-D reconstructions were performed on a tetrahedral mesh discretizing the mouse model with 7712 nodes.

Fig. 3.

Forward results of spherical targets with different radii: (a) 1.0 mm, (b) 1.5 mm, (c) 2.0 mm, and (d) 2.5 mm.

3.3 Reconstruction Results

Based on the joint $\ell_{1}$ and Laplacian manifold regularization model and the experiment settings, we utilize the GPSR method to solve the inverse problem. Enough iterations were used to ensure sufficient convergence. The regularization parameters of the method were manually optimized to balance the reconstruction error and image contrast.

Specifically, parameter $\tau$ and $\lambda$ for Laplacian regularization term have significant influence on the reconstructed results. This is mainly because the variable $X$ and input data $\Phi$ are $N \times 4$ matrices, for each column the optimal parameters are different. Thus, it is necessary to add constraints to the parameters to get the best results of each column at once.

Figures 4 and 5 show the reconstruction results in 3-D views and cross-sectional views respectively. It is clearly seen from the reconstructed images at different time points that the center of reconstructed fluorophores is close to the spherical targets and the quality is acceptable. The quantitative assessment of the reconstruction results are shown in Table 2.

Fig. 4.

Reconstruction results. (a), (b), (c) and (d) shows the reconstruction results of spherical targets with radius equals to 1.0, 1.5, 2.0, 2.5 mm, respectively.

Fig. 5.

The cross-sectional views of the reconstruction results. (a)–(d) are correspond to Fig. 4. (a)–(d). The black circle is the real position of target.

Table 2. Quantitative assessment

Table 2 shows that the distance between the center of the reconstructed fluorophore and real target is 1.1 mm to 1.5 mm, which proves the reliability of the proposed method. Besides, in consideration of the radius of target varies from 1.0 mm to 2.5 mm, the proposed algorithms shows stability in reconstructing small target and big target.

Figure 6 shows the volume of spherical targets and reconstructed targets. The transvers axis represents time and the vertical axis represents value of volume. With the increase of radius, the reconstruction results show growth trend, and the differences between the real targets and reconstructed targets are relatively small.

Fig. 6.

The volume of spherical target and reconstructed target. The curves are depicted to different colors (spherical target: dark yellow; reconstructed target: grey) (color figure online).

4 Discussion and Conclusion

In this paper, we propose a joint $\ell_{1}$ and Laplacian manifold regularization model for dynamic FMT reconstruction. Compared with other methods, the proposed approach makes full use of the sparsity and spatial structure information of the boundary measurement data, and a sequence of FMT images of different time points can be reconstructed in one step. Besides, numerical simulation is implemented to validate the performance of the joint regularization model and the proposed algorithms. The results demonstrate the reliability and stability of the proposed approach. In general, this approach for dynamic FMT can provide satisfactory reconstructed images.

In the simulation section, we mainly focus on testing the proposed model and the algorithm, only the first experiment, assuming that all the fluorescent targets have the same concentration, was illustrated. The other comprehensive experiments will be carried out in future works.