Shape-based reconstruction of dynamic fluorescent yield with a level set method

Background Fluorescence molecular tomography (FMT) is an optical imaging technique that reveals biological processes within small animals through non-invasively reconstructing the distributions of fluorescent agents. The primary problem in FMT with non-stationary fluorescent yield is the increase of the unknown parameters to be reconstructed. In this paper, a method is proposed to reconstruct dynamic fluorescent yield. Methods A shape-based reconstruction method that recovers dynamic fluorescent yield with a level set method is proposed for FMT. To reduce the number of unknown parameters, a level set function is introduced to describe the shape of target and a small number of parameters are used to describe the fluorescent yields at different time points. Results Results of simulations and phantom experiments demonstrate that the proposed method can recover well the dynamic fluorescent yields, shapes and locations of the target. Conclusions The proposed method can handle the cases with non-stationary fluorescent yields and recover the fluorescent yields at each projection angle.

reconstruct a static distribution of fluorescent yield from data sets acquired within several minutes. The other kind of methods is based on different models. These methods model the metabolic process with a two-compartmental model [7,8] or a polygonal line model [9,10], and reconstruct the distributions of four pharmacokinetic parameters in the two-compartmental model or the distributions of fluorophore concentration variation rate, instead of the distribution of fluorescent yield. These model-based methods convert the inverse problem of dynamic FMT from reconstructing a series of static distributions of fluorescent yield to reconstructing the distributions of other parameters that can represent the metabolic information. Their limitations are the restriction of the used model. In this paper, a new approach is proposed to reconstruct dynamic fluorescent yield based on a level set method. Different from the model-based methods, this method still reconstructs the distributions of fluorescent yield, thus the target metabolic curves are not restricted by specific models. This method is capable of recovering the distributions of fluorescent yield at each projection angle.
If the fluorescent yield is considered as the unknown parameter, the primary problem in FMT reconstruction with dynamic fluorescent yield is that the number of unknown parameters to be reconstructed is dramatically increased, because the fluorescent yield varies projection by projection in the dynamic problem whereas it does not change in the static problem. If the distributions of fluorescent yield at each projection angle are considered to be reconstructed, the number of unknown parameters will be multiplied by the number of projection angles (N p ). Based on prior information that the same tissue or organ should have similar metabolic property, a shape-based reconstruction strategy is implemented in FMT to reduce the number of unknown parameters. In this strategy, the reconstructed object is divided into numbers of disjoint subregions and the shapes of these subregions are represented by a level set method [11] or spherical harmonics expansion [12]. Accordingly, the fluorescent yield in each subregion is described by only one parameter. It changes the reconstruction problem from reconstruction of fluorescent yield at each voxel or node to reconstruction of level set function or expansion coefficients with several parameters of fluorescent yield. When a problem with dynamic fluorescent yield is considered, the number of unknown parameters will be reduced because the shapes of the distributions of fluorescent yield at different projection angles are described by the same set of unknown parameters by using the level set method or spherical harmonics expansion. In this paper, the level set method is chosen to represent the shapes of subregions for its ability of arbitrary boundary representation. Although spherical harmonics expansion can also describe arbitrary boundary theoretically, the maximum degree of spherical harmonics, which controls the shape representation, needs to be determined manually [13].

Methods
The generation and propagation of fluorescence in continuous wave FMT can be described by a couple of diffusion equations with Robin-type boundary condition, which are given by [14,15]: where Φ is the photon density. r denotes the coordinates and Ω denotes the region of imaged object. The subscripts x and m represent excitation and emission, respectively. Q is the excitation source term. n is the outward normal vector to the surface and q is a term related to the optical reflective index mismatch at the boundary. μ a and D are the absorption and diffusion coefficients, respectively. η is the quantum yield and μ af is the absorption coefficient of fluorophore to the excitation light. Their product is the fluorescent yield which needs to be reconstructed. Numerical approaches are commonly used to solve the diffusion equations with arbitrary shaped domains. In this work, finite element method (FEM) is chosen to solve the forward problem, which converts the coupled diffusion equations given by Eq. (1) into the following linear equations after discretizing the imaged domain into a mesh with N n nodes [16]: Here Φ, X, and Q are column vectors with a length of N n • X = ημ af is the vector of the unknown fluorescent yield. K is an N n × N n matrix which is commonly called stiffness matrix in FEM. F is a matrix with the same size as K, which is related to the excitation light field. K and F are constructed according to the topological structure of the mesh. K, F, and Q are given by: where the subscripts i and j denote the indices of column and row, respectively. φ i is the shape function at node i.
To reduce the number of unknown parameters, the shape-based reconstruction strategy is used. The imaged domain is assumed to consist of several disjoint subregions and the fluorescent yield in each subregion is represented by a single parameter. For simplicity, the case with two subregions (target region R f and background region R b ) is considered here and x f and x b are used to denote the corresponding fluorescent yields. To describe the shape of the subregions, a level set function is introduced [17]: where ψ is the level set function. Because the case with dynamic fluorescent yield is considered, x, x f , and x b are functions of time. In order to combine the level set function with FEM, ψ is discretized into an N n dimensional vector Ψ of basis coefficients by a basis expansion with the shape functions as follows: where ψ i is the level set function at node i.
Then Eq. (7) is reformed by a single equation as follows: where H (ψ) is the step function which is equal to 0 or 1 if ψ ≤ 0 or ψ > 0. Through Eqs. (8) and (9), the unknown parameters to be reconstructed are converted from numbers of vectors X into only one vector Ψ with a set of values of fluorescent yields {x f (i)} and {x b (i)} after the discretization of the time t.
To reconstruct Ψ, {x f (i)} and {x b (i)} from boundary measurements, the following object function is defined: where P is the measurement operation vector [18] used to extract the photon density at a certain detector point from the vector Φ. N f denotes the number of frames of data used in the reconstruction and N d denotes the number of detector points at each projection angle. y * and y are the predicted measurements calculated from the forward problem and measurements acquired by the system, respectively.
Differentiating Eq. (10) with respect to ψ m yields:   where the superscripts (n) and (n + 1) denote the indices of iteration. τ 1 and τ 2 are the step lengths. They are determined empirically according to the value of the cost function. Y k, i is the vector of measurements acquired at the kth frame and ith projection angle. J is the Jacobian matrix of the forward operator, which consists of row vectors {P ij K −1 F i }. Because the Dirac function in Eq. (11) cannot be numerically calculated, it is omitted in the reconstruction. During the reconstruction, a low-pass filter [20] is applied to the reconstructed fluorescent yields after each iteration to smooth the curve of the reconstructed fluorescent yields.

Results and discussion
Simulation and phantom studies were implemented to validate the performance of the proposed method. The experimental setup of the simulation and phantom studies is shown in Fig. 1. The imaged object was a cylinder with a diameter of 3 cm and height of 5.2 cm as shown in Fig. 1a. A tube with a diameter of 0.4 cm and height of 5.2 cm was inserted in the cylinder as the target. The tube was located at the position (x, y) = (0, 0.5) cm. In the simulations, a varied fluorescent yield according to the red  Fig. 1b was assigned to the tube. The curve was generated with a metabolic model of indocyanine green (ICG) [21]. Measurements were generated by Eq. (3) and 1 % Gaussian noise was added. In the phantom experiments, ten sets of measurements were acquired when different concentrations of ICG according to the blue points in Fig. 1b were filled in the tube. Then the dynamic measurements were synthesized by linear interpolation. Simulations and phantom experiments on a model with double targets as shown in Fig. 2a were also implemented. Two tubes were inserted in the cylinder as the fluorescent targets and the fluorescent yields or concentrations of ICG were set according to two different curves as shown in Fig. 2b. The heights and radii of the cylinder and tubes were the same as those in the previous experimental setting. The two tubes were away from each other with an edge-to-edge distance of 0.8 cm.
In the phantom experiments, 1 % intralipid with μ s ′ = 10.0 cm −1 and μ a = 0.02 cm −1 was filled in the transparent cylinder. Accordingly, the same optical coefficients were used in the simulations. A parallel excitation based FMT system [22] as shown in Fig. 3 was used to acquire datasets in the experiments. The system consisted of three parts: excitation module, rotation stage, and detection module. The excitation module was a 300 W xenon lamp (Asahi Spectra, Torrance, CA, USA) coupled with a 770 ± 6 nm bandpass filter. A slit was used to generate a line illumination as the excitation source. The detection module consisted of a −70º cooled electron multiplying charge-coupled device (EMCCD) camera (iXon DU-897, Andor Technologies, Belfast, Northern Ireland, U. K.), a Nikkor 60 mm, f/2.8D lens (Nikon, Melville, NY, USA), and an 840 ± 6 nm bandpass filter. The phantom was placed on the rotation stage between the excitation and detection modules and measurements were acquired at 18 projection angles over 360º illuminations and detections. The time for the data acquisition of each projection angle was 2 s. The same time interval was used in the simulations. After the detection of fluorescence, 72 white light images were acquired to recover the geometry of the imaged object with a surface reconstruction method [23]. In the reconstructions, the initial level set functions Ψ (0) were set as a small positive constant for every node, but the value (i.e., 1) used in the phantom studies was chosen to be larger than the one (i.e., 0.001) used in the simulations because the data acquired by EMCCD and the data used in the simulations had different orders of magnitude. Since the level set functions were set to be homogenous at the beginning, the initial subregions consisted of only the background. The initial values of the fluorescent yield x b were set as 0 for all time points. It should be noticed that x f and x b cannot be the same, otherwise the gradient given by Eq. The time interval between two neighbor fluorescent yields is equal to the time for the data acquisition of each projection angle (2 s). It can be found that the reconstructions of both simulations and phantom experiments recover well the dynamic fluorescent yields, shapes and locations of the target. The reconstructed time course for the phantom experiments shows greater mismatch with the true one than that for the simulations, which is probably due to the complicated noise in the experimental datasets of different frames.
In this paper, the iteration formulas given by Eqs. (13) and (14) are derived from the gradient descent method. These formulas are similar to those derived from an artificial time evolution approach [11,17,24]. The Dirac function in the gradient for the level set function is omitted in both the proposed method and the methods based on the artificial time evolution approach. Rigorously, Eq. (13) is only proper for those locations with ψ = 0, because δ (ψ) is not zero only when ψ = 0. However, the gradient given by Eq. (13) shows the capability of finding the subregion of the target by decreasing the level set function with positive values until it becomes negative in numerical applications.

Conclusion
In general, a shape-based reconstruction method that recovers dynamic fluorescent yield with a level set method is proposed in this paper. This method is superior to the conventional method that reconstructs only a static distribution of fluorescent yield for the reconstruction of dynamic FMT, because it is capable of resolving the fluorescent yields