A Sparsity-Constrained Preconditioned Kaczmarz Reconstruction Method for Fluorescence Molecular Tomography

Fluorescence molecular tomography (FMT) is an imaging technique that can localize and quantify fluorescent markers to resolve biological processes at molecular and cellular levels. Owing to a limited number of measurements and a large number of unknowns as well as the diffusive transport of photons in biological tissues, the inverse problem in FMT is usually highly ill-posed. In this work, a sparsity-constrained preconditioned Kaczmarz (SCP-Kaczmarz) method is proposed to reconstruct the fluorescent target for FMT. The SCP-Kaczmarz method uses the preconditioning strategy to minimize the correlation between the rows of the forward matrix and constrains the Kaczmarz iteration results to be sparse. Numerical simulation and phantom and in vivo experiments were performed to test the efficiency of the proposed method. The results demonstrate that both the convergence and accuracy of the proposed method are improved compared with the classical memory-efficient low-cost Kaczmarz method.


Introduction
Fluorescence molecular tomography (FMT) is an imaging modality that can localize and quantify fluorescent markers to resolve biological processes at molecular and cellular levels. Being minimally invasive, relatively inexpensive, and portable, FMT has been successfully applied in small animal research and preclinical diagnostics such as cancer diagnosis, drug discovery, and gene expression visualization [1][2][3][4].
Due to a large number of unknowns and a limited number of measurements as well as the diffusive transport of photons in biological tissues, FMT reconstruction is an ill-posed inverse problem [5][6][7]. To improve the FMT imaging quality, both the noncontact FMT technique [8,9] and the strategy of multiple excitations can be used to obtain more measurements. Structural a priori information provided by CT or MRI can also be incorporated into FMT imaging [3,[10][11][12]. Moreover, reconstruction algorithms can resort to regularization strategies and find meaningful and numerically stable solutions. In [13,14], the Tikhonov regularization, namely, presented in [18]. In [19], a hybrid regularization method incorporating 1 and 2 norm penalty is proposed to recover the 3D fluorophore distribution. In these techniques, optimal selection of the regularization parameter is needed to avoid over-or underregularization. Being a memory-efficient lowcost numerical solver that avoids bulky matrix computations in large-scale problems, Kaczmarz algorithm, also known as algebraic reconstruction technique (ART), iteratively updates the solution using only one equation at a time and has been applied in optical tomographic reconstruction [20][21][22]. During reconstruction, the Kaczmarz method may use the measurements in the order that they are collected, which is known as the sequential access order. To speed up the convergence rate of the Kaczmarz method and give better results in the first iteration relative to the sequential access scheme, different access orders have been proposed [23][24][25]. The idea of these different access orders is to minimize the correlation between measurements that are successively accessed by the iterative projection inversion method. In [20], the influence of the data access order is investigated when Kaczmarz method is used to perform diffuse optical tomography. The study shows that the convergence speed can be significantly improved by selecting proper projection access order.
In FMT, the forward matrix maps the fluorescent targets to the surface measurements. Generally, the rows of the forward matrix are correlated because of the correlations among source-detector maps from the same projection and the interrelations of different projections [26,27]. In this work, we present a strategy which computes a preconditioning matrix to minimize the coherence of the preconditioned forward matrix. Then the Kaczmarz method which uses the sequential access order is adopted to solve the preconditioned FMT reconstruction problem. After preconditioning, the projections are close to perpendicularity and the convergence rate of the Kaczmarz method can be speeded up. As most optical fluorophores are designed to accumulate in relatively small, specific regions in tissues, such as tumors, and hence the fluorophore distributes sparsely in the imaging domain, we propose sparsity-constrained reconstruction method to perform FMT and the method is named as sparsity-constrained preconditioned Kaczmarz (SCP-Kaczmarz) method. Different from the existed 1 norm regularization methods, this proposed SCP-Kaczmarz method adopts a thresholding step to the Kaczmarz results to satisfy a user-defined sparsity value.
The remaining of this paper is organized as follows. We first describe the mathematical forward model for FMT imaging, then the SCP-Kaczmarz method is presented for FMT reconstruction, then the numerical simulation and physical phantom and in vivo experiments are performed to evaluate the proposed method, and finally the discussion and conclusion are given.

Forward
Model for FMT Imaging. When a CW point laser is used as excitation light, the diffusion of excitation and emission lights through biological tissues can be described by two coupled diffusion equations with the Robin-type boundary condition, and the coupled diffusion equations can be presented as follows [28]: where ∈ Ω, Ω being the domain under consideration. The subscripts and denote excitation light and emission light, respectively. , = 1/3( , + (1 − ) , ) is the diffusion coefficient with , as the scattering coefficient, is the anisotropy parameter, and , is the absorption coefficient. Φ , denotes the photon density. The fluorescent yield is the unknown parameter to be reconstructed, which is denoted as hereafter. By using the finite element method (FEM), the linear relationship between the boundary measurements and the desired unknown fluorophore distribution can be obtained from (1) and is described by where A is the forward matrix, the sizes of , A, and are × 1, × , and × 1, respectively. is the number of surface measurements and is the number of unknowns needed to be determined inside the imaging domain. Usually < , and this means that the number of measurements is smaller than that of the unknowns.

Sparsity-Constrained Preconditioned Kaczmarz Method.
It is known to us that the convergence of the Kaczmarz method is affected by the data access order. If the measurements are prearranged in a scheme that the projections are close to perpendicularity, the convergence of the Kaczmarz method will be speeded up. In this paper, rather than changing the sequential data access order, we design a preconditioner to minimize the correlation between the rows of the forward matrix of FMT and hence to make the Kaczmarz method converge quickly. Denote the preconditioning matrix as W and the preconditioned forward matrix as B, then we get B = WA, and W can be obtained by solving the following optimization problem: where I is the × identity matrix and ‖ ⋅ ‖ F is the Frobenius norm.
Considering the singular value decomposition of A which is described by Equation (4) means that the rows of the preconditioned forward matrix are orthogonal to one another. If the preconditioner W is badly conditioned, we can use the diagonal loading strategy to mitigate the ill condition and W is calculated by W = (SS T + I) −1/2 , where ≪ 1 is a constant [27]. Figure 1 provides a geometric insight into the iterative progress of the Kaczmarz and the preconditioned Kaczmarz algorithms. Figure 1(a) presents a geometrical interpretation of Kaczmarz applied to a 2D problem. Here, each line represents a hyperplane in the solution space corresponding to one of the equations, and the solution is the intersection of the dashed lines. The progress of Kaczmarz is represented by dark blue dots and arrow lines. As depicted in Figure 1(a), the points with dots iteratively progress toward the solution (intersection of the two dashed lines) by orthogonal successive projections onto the two lines [22]. In Figure 1(b), the blue diamond and arrow demonstrate the convergence of preconditioned Kaczmarz algorithm toward the solution. Because the forward matrix has been preconditioned, the two green dashed lines which demonstrate the hyperplanes corresponding to the two equations are perpendicular. In theory, only one iteration is needed for the algorithm to converge to the solution. However, because of the ill condition of the forward matrix in FMT imaging and the presence of noise, the two lines are not completely perpendicular. Multiplying (2) by W on both sides, we can obtain the distribution of the fluorescent targets by solving We use the classical Kaczmarz technique to solve (5) and the unknown is updated by where B is the th row of B.
Considering that the fluorescent target is sparsely distributed, we add a thresholding step to (6) to make the sparsity of the Kaczmarz result remain as close as possible to a preset value. The thresholding procedure is depicted by where ( ) is the th element of ( ) and is between 0 and 1 and can be obtained by solving the minimization problem [29] = arg min sparsity ( ( ) ( )) − .
The thresholding step also guarantees the nonnegativity of the solution. The minimization problem of (8) can be solved by using linear searching strategy. And the sparsity of in (8) is defined as [30] sparsity where is the size of vector and ‖ ‖ 1 and ‖ ‖ 2 denote 1 and 2 norm of , respectively. The curve in Figure 2(a) depicts the variation of sparsity value with the number of nonzero elements in (assume that has 300 elements and the nonzeros in are constant, e.g., 1), from which we can see that the sparsity value ranges from 0 for nonsparse results to 1 for extremely sparse results. Figure 2(b) shows with sparsity value of 1, 0.87, and 0.73, respectively. So far, the implementation of the proposed SCP-Kaczmarz method can be summarized as follows: (1) Initialize ( ) , where = 0, and preset the wanted sparsity value . (2) Perform singular value decomposition to A and compute the preconditioning matrix W. (4) Update ( ) from ( −1) by solving = B using the classical Kaczmarz method. (5) Keep the large elements of ( ) and set the other elements to zero to make the sparsity of the result equal to the wanted sparsity value .
A parameter, named as the wanted sparsity value , should be predetermined for the proposed SCP-Kaczmarz method. As we do not know the true distribution of the fluorescent target in practice, the ratio of the volume of the fluorescent target to that of the imaging domain can be first estimated, and then the corresponding sparsity is calculated as the wanted according to (9) under the assumption that the fluorescent target is uniformly distributed.

Experiments and Results.
To demonstrate the performance of the proposed method, numerical simulation and phantom and in vivo experiments were conducted. In these experiments, cases with single fluorescent target and multiple fluorescent targets were considered. All the reconstructions were implemented on a personal computer with an 8 GB memory and an Intel-Core i7 CPU. The relative deviation, Dice coefficient, and sparsity value were calculated to evaluate the SCP-Kaczmarz algorithm. The relative deviation is defined by = ‖ − ‖ 2 /‖ ‖ 2 , where is the reconstructed target and is the actual target. The Dice coefficient is where ⋅ is Hadamard product and the sparsity is defined by (9). In addition, the computational time and memory consumption of both the methods were also recorded.

Numerical Simulation Experiments.
In the numerical simulation, a 3D digital mouse atlas of CT and cryosection data was utilized to provide anatomical information [31]. Cases with single, two, and three fluorescent targets inside the mouse atlas were studied, respectively.

Reconstruction of Single Fluorescent Target.
In this section, one fluorescent target inside the digital mouse atlas was reconstructed and two imaging models were investigated. The first model is reconstruction of small fluorescent target while the second one is reconstruction of big fluorescent target. The small target model is usually used to mimic small tissue with fluorescent probe, such as the early tumor; the big target model can be used to recover the biodistribution of fluorescent probe in organs, which is important in drug pharmacokinetics study [32]. The small target model is shown in Figure 3(a), where a spherical fluorescent target (marked in red color) with radius of 1.5 mm was placed at "12.9 mm, 9.9 mm, and 16.5 mm" in the liver. The optical parameters of the mouse organs including muscle, heart, lungs, liver, kidneys, and stomach were listed in Table 1 [19]. As illustrated in Figure 3(b), the fluorescent target was excited sequentially by 5-isotropic point sources located one mean free path of photon transport beneath the mouse surface on = 16.5 mm plane. For inverse reconstruction, the atlas torso was discretized into 38735 tetrahedrons and 7511 nodes. The sparsity value and fluorophore distribution were set to be 0.9 and 0, respectively, when = 0. Figure 4 shows the relative deviation, the Dice coefficient, and the sparsity value obtained during the iteration process of the SCP-Kaczmarz and Kaczmarz algorithm. From Figure 4 we can see that, after about 20 iterations, the proposed method converges to the true solution with sparsity value of 0.9476. Figure 5 shows the 3D targets and 2D slices at = 16.5 mm recovered by the SCP-Kaczmarz and Kaczmarz method after 100 iterations, respectively. The actual target is indicated by the black circle and the reconstructed results are normalized by the true intensity. It can be seen that both methods can locate the target accurately, while the first one is more appropriate for quantitative analysis and profile reconstruction. It is known to us that the singular value decomposition of the forward matrix, needed for the proposed method to calculate the preconditioner, is computationally expensive. Fortunately, the SVD only needs to be performed for one time and can be done before the iteration starts. So, in this work, we just measure the elapsed time for the proposed algorithm to iteratively solve the preconditioned FMT inverse problem by using the MATLAB functions tic and toc and the SVD time cost is not included. The time cost is 79 seconds and 100 seconds, respectively, for the SCP-Kaczmarz method and the Kaczmarz method to run 100 iterations. The corresponding memory consumption, which is calculated by using the MATLAB instruction profile on memory, is 565040 KB and 565016 KB. As the memory is mainly used to store the preconditioned or the original forward matrix, the two algorithms have similar memory cost.
The big target model was used to recover the biodistribution of fluorescent probe in heart. The target was excited sequentially by 5-isotropic point sources located at one mean free path of photon transport beneath the mouse surface on the = 7.3 mm plane where the heart centered. Figure 6 illustrates the reconstruction results (which are normalized to the intensity of the actual target) obtained by the SCP-Kaczmarz and the Kaczmarz method after 100 iterations, where the actual heart inside the body is hidden for clarity. We also plotted the recovered intensity at each node of the discretized atlas torso, as shown in Figure 7, where theaxis denotes the node index and the -axis denotes the intensity. The relative deviation, Dice coefficient, and sparsity are illustrated in Figure 8. Although the two methods get different results, the sparsity values are the same after 40 iterations. The computational time of the SCP-Kaczmarz and Kaczmarz is 80 seconds and 81 seconds, respectively, and both the memory usages are about 565 MB.

Reconstruction of Multiple Fluorescent Targets.
In this section, we used the proposed SCP-Kaczmarz method to recover multiple targets. As shown in Figure 9(a), two spherical fluorescent targets with radius of 1 mm were placed at "13 mm, 12 mm, and 16.5 mm" and "13 mm, 6 mm, and 16.5 mm" in the liver of the digital mouse. Figure 9 To further test the ability of the proposed method to distinguish multiple targets, we considered three fluorescent spheres with radius of 1 mm placed in the liver of the digital mouse. As shown in Figure 10

Phantom and In Vivo Experiments.
In this section, our homemade dual-modality FMT-Micro-CT imaging system [19] was used to perform the phantom and in vivo experiments. Two phantom experiments were conducted. In the first experiment, a 20 mm cubic phantom made from polyoxymethylene was placed on the rotational stage of the imaging system. A small hole with 1 mm radius and 2 mm height was drilled at "15 mm, 7 mm, and 9.  = 0.0114 mm −1 and = 10.1 mm −1 , respectively. In the second experiment, a cylinder phantom with 10 mm radius and 30 mm height was used. Two holes with 1 mm radius and 5 mm height were drilled in the cylinder phantom. The two holes were centered at "5 mm, 4 mm, and 15 mm" and "5 mm, −4 mm, and 15 mm," respectively. Both holes were filled with 40 M of Cy5.5 solution. The 671 nm CW laser was used as the point source. Five excitation points were set equivalently along the right half side of the phantom on = 15 mm plane. The cubic phantom is illustrated in Figure 11(a), where the green cylinder denotes the Cy5.5 solution. The initial sparsity value was set to be 0.8 for the first iteration. Figures  11(b) and 11(c) show the normalized results at = 9.5 mm by the SC-Kaczmarz and the Kaczmarz method after 200 iterations, respectively, where the true target is indicated by the white circle and the location errors are about 0.7 mm and 1.9 mm, respectively. Compared with the Kaczmarz result, the intensity of the SCP-Kaczmarz result is larger due to the fact that the latter one distributes more sparsely. The time cost is 55 seconds and 58 seconds and the memory cost is about 100 MB.
The cylinder phantom is shown in Figure 12 for the SCP-Kaczmarz algorithm to iterate 100 times and 862 seconds is needed for the Kaczmarz algorithm to iterate 500 times. The memory used by the two methods is about 600 MB, where the phantom is discretized into 7851 nodes for inverse reconstruction.
The proposed SCP-Kaczmarz method was also used to recover the fluorescent target from the in vivo small animal experimental data [32]. In the experiment, the fluorescent target was made of a glass tube full of 4000 nM Cy5.5 solution and was implanted into the abdomen of an adult BALB/C mouse. A 671 nm CW point laser was used to excite the target at four positions sequentially and the optical signal on the mouse surface was collected. After the optical images acquisition, the anesthetized mouse was scanned by the Micro-CT subsystem. The reconstructed CT slices were segmented into five components (heart, lungs, liver, kidneys, and muscle) and used to provide prior structural information for the FMT inverse problem. The sparsity value was set to be 0.7 for the first iteration. The in vivo experiment results are shown in Figure 13

Discussion and Conclusion
In this work, we propose a sparsity-constrained preconditioned Kaczmarz method to solve the inverse problem in FMT. First, a preconditioner is computed to minimize the correlation between the rows of the FMT forward matrix, then the classical Kaczmarz method is used to solve the preconditioned inverse problem, and finally the large elements of the Kaczmarz solution are kept and the other elements are set to zero to make the result satisfy a preset sparsity value. The threshold value is obtained by solving a minimization problem using linear searching strategy. The performance of the proposed algorithm is demonstrated by numerical simulation and phantom and in vivo experiments. In numerical simulation, both small target and big target can be recovered with high accuracy by the proposed method. As the correlation between the rows of the forward matrix has been minimized, the SCP-Kaczmarz method converges to the true solution rather faster than the classical Kaczmarz method. In the phantom and in vivo experiment, the proposed algorithm shows performance improvement both in location accuracy and in convergence speed relative to the classical Kaczmarz technique. Furthermore, experiments with two targets and three targets were performed. The results show that both the methods can get satisfactory results. However, we find that the Kaczmarz method converges rather slowly when recovering more than one target, while the convergence of the proposed SCP-Kaczmarz is not affected by the number of targets. The reason may be that the correlation between the rows of the forward matrix when multiple targets are present is stronger compared to the case when one target is present. And the stronger correlation leads to slower convergence of the Kaczmarz method. After preconditioning, the correlation is minimized and hence the convergence remains the same for different number of targets. The computational time and the memory usage are also calculated. The results show that, under the same number of iterations, the Kaczmarz method runs a little faster than the SCP-Kaczmarz method, and they consume similar memory. That is because the memory is mainly used to store the preconditioned or the original forward matrix which are of the same size. A parameter, named the wanted sparsity value, should be predetermined when using the SCP-Kaczmarz method. As the number of nodes in the imaging domain is known, we can estimate the number of nodes the fluorescent target covers and hence get a sparsity value by (9). However, we found that the sparsity value has little effect on the reconstruction result in our experiment. It is known to us that the singular value decomposition of the forward matrix, needed for the proposed method to calculate the preconditioner, is computationally expensive. Fortunately, the SVD only needs to be performed for one time and can be done before the iteration starts. No additional computation load will be caused when the proposed method iteratively solves the preconditioned FMT inverse problem. In numerical simulation, as the surface data was obtained based on the diffusion approximation (DA) model, which is also used in solving the inverse problem, the proposed algorithm performs well to reconstruct the intensity and the shape of the fluorescent target. For real experiment, we do not know actually the light transport model. So there is a mismatch between the true model and the diffusion approximation, and the proposed method does not perform as well as it does in the simulation experiment. In the future, we will focus on light transport model based on higher order approximation (e.g., the SP approximation) and perform in vivo experiment with multiple targets to further investigate the proposed method.