A Fast Reconstruction Algorithm for Fluorescence Optical Diffusion Tomography Based on Preiteration

Fluorescence optical diffusion tomography in the near-infrared (NIR) bandwidth is considered to be one of the most promising ways for noninvasive molecular-based imaging. Many reconstructive approaches to it utilize iterative methods for data inversion. However, they are time-consuming and they are far from meeting the real-time imaging demands. In this work, a fast preiteration algorithm based on the generalized inverse matrix is proposed. This method needs only one step of matrix-vector multiplication online, by pushing the iteration process to be executed offline. In the preiteration process, the second-order iterative format is employed to exponentially accelerate the convergence. Simulations based on an analytical diffusion model show that the distribution of fluorescent yield can be well estimated by this algorithm and the reconstructed speed is remarkably increased.


INTRODUCTION
With the discovery of biocompatible, specific fluorescent probes and the development of imaging technologies, the potential of fluorescence tomography as a means for molecularly based noninvasive imaging of biological tissues has received in recent years increased attention [1][2][3]. Fluorescent beacons emitting in the near-infrared (NIR) bandwidth are always preferred, since hemoglobin and water absorb minimally in this spectral window so as to allow photons to penetrate for several centimeters in tissues [4].
Using preferentially accumulated fluorescent probes as indicators or contrast agents, fluorescence optical diffusion tomography (FODT) is performed by launching light at the probes' excitation wavelength into the tissue. The fluorescent beacon absorbs the incident light, and emits light at a longer wavelength when it drops to the ground state. Then the emission is measured by an array of detection devices at the surface of the body. However, as the strong diffusion of NIR in biological tissues, reconstruction of very large unknown inside characteristics from the limited detected data at the boundary is one of the main difficulties in FODT. Many reconstructive approaches utilize iterative methods for data inversion, such as the algebraic reconstruction technique (ART) [5], Newton's or Newton-type optimization methods [6,7], and Bayesian nonlinear least-square method [8,9]. They are always time-consuming and far from meeting the real-time imaging demands.
In this study, a fast algorithm based on the preiteration is applied to the inversion process of fluorescence tomography. For simulating the photon's propagation in tissues with fluorescent beacons inside, a previously reported DPDW model based on Born approximation is simply introduced at the beginning. Then, the preiteration fast algorithms are presented in detail, emphasizing the second-order method. After that, the simulation using the second-order form is investigated and the results are shown. Finally, we analyze the computation burden and convergence property of the second-order iteration form and give the conclusion.

DPDW MODEL
Often a couple of diffusion equations in frequency-domain is employed to describe the propagation of both excited light and fluorescent light in diffusive medium, that is [6,7,10] where Φ x,m is the photon density for excitation (subscript x) or fluorescent light (subscript m), D x,m (r) is the diffusion coefficient, and μ ax,m (r) is the absorption coefficient. Based on this model, the fluorescence lifetime τ(r) and the yield η(r) can be estimated through the boundary measurements.
Equation (1) can be solved by analytical or numerical methods. In this paper, we use an analytical model of Born approximation for specific medium geometry to demonstrate the inversion algorithm. In fact, the fast algorithm could also be applied to arbitrary geometries, where the model is discretized by numerical methods [7,10] or the Kirchhoff approximation [5].
In the frequency-domain model, an amplitude-modulated incident point source of photons into a diffusive medium produces a diffuse photon density wave (DPDW) [11][12][13]. Let an intensity-modulated point source of amplitude Θ 0 be located at r s in a homogeneous infinite medium. Then the spatial part of the originating DPDW at position r is [11] with the wave number k = [(−νμ a + iω)/D] 1/2 , and D = ν/3μ s is the diffusion coefficient with the reduced scattering coefficient μ s and the speed of light in the medium ν. Here, ω is the angular modulation frequency of the source. Treating fluorescent beacons as two-level quantum systems and assuming that there are no saturation or photon quenching effects, the fluorescent photon density δu f l , measured at a detector position r di due to a localized probe with volume d 3 r k embedded within the medium, is [11] δu f l r k , r s j , r di with the excited source at r s j . Here, λ 1 and λ 2 represent the excited light wavelength and the fluorescent wavelength in the near-infrared section, respectively. G(r d − r, k λ2 ) = exp(ik λ2 |r d − r|)/4π|r d − r| is Green's function solution to the diffusion equation and represents the variance of fluorescent DPDW from fluorescent probe to the detector. For a weakly absorbing spatial distribution of fluorescent probes, the detected fluorescent DPDW at r di can be found by integrating overall fluorescent sources [11,12]. In the reconstruction, for the measurement at positions r di (i = 1, 2, . . . ,M i ), the integral can be digitized as due to the sources r s j ( j = 1, 2, . . . , M j ). As only one of the sources is working at a time, the total number of mea- In fluorescence tomography, continuous-wave (CW) mode is always chosen, that is, ω = 0, and only η is reconstructed. Then substituting (2) in (3) will lead to the following matrix equation: where U represents an M × 1 column vector of the detected data, X is a column vector of unknown values of fluorescent yield η at N reconstructed points, and matrix A indicates the obtained M × N weighted coefficients.

PREITERATION INVERSE ALGORITHM
As in FODT, the inside reconstructed points number N is always much bigger than M, the measurement number at the boundary, the equation series (4) is always ill-posed and indefinite. In this case, the direct inverse matrix of A does not exist. However, its generalized inverse can be employed to solve (4).

Preiteration algorithm based on generalized inverse
If the Moore-Penrose inverse of A exists and is known as A + , the unique solution of (4) which has the minimum norm and the least square can be obtained simply by [14] There are several direct methods to calculate the generalized inverse A + , for example, regularized SVD method. However, the iterative method is always preferred in computerized calculation, especially for large datasets, as it is easy to be programed and occupies much less ram than direct methods. Supposing the residual error series R k = I − A S k (I is the unit matrix of M × M), series will be convergent to A + when k → ∞ [14]. Here S 0 can be chosen as αA T [15], with α = 1/λ max . And λ max is the maximum eigenvalue of A·A T , where A T is the transposed matrix of A. From the analysis above, a two-step reconstructed algorithm can be formed.
(1) Offline preiterative step: the approximation of generalized inverse A + is calculated by several iterative steps of (6). (2) Online reconstruction: when the weighted matrix A keeps unchanged or the variation can be ignored, for updated detection U the unknown character X can be reconstructed simply through (5).
This preiteration method has already been applied to the image reconstruction in electrical impedance tomography (EIT) [15] which also belongs to the so-called "soft field" imaging as FODT, and it is proved that Landweber iteration method, which can produce higher quality reconstructed image than other direct regularized methods, is in fact a modification of the above preiteration algorithm [16]. However, Xiaolei Song et al. 3 compared with Landweber method, the preiteration method remarkably improves the reconstructing speed by performing the time-consuming iterative process offline.

Second-order iteration form
However, the first-order preiteration algorithm with form equation (6) needs the same iteration steps as the Landweber method to produce the same quality images [15]. So just like the slow convergence of Landweber, for larger-sized dataset in FODT, iteration form of (6) is also very time-consuming even in the preiteration process. In order to speed up the iteration process, the second-order iterative format is used in our work.
To prove the convergence of the second-order form equation (7), we examined the convergence of S k and the residual error R k as follows.
First, by including (7), the iterative formula of R k can be obtained as Then it can be inferred that According to (7) and (9), S k+1 can be written as a function of R 0 and S 0 in the formula if ρ(R 0 ) < 1, let k → ∞, it yields where S ∞ can be proved as the generalized inverse matrix of A [14]. However, with the first-order iteration form equation (6), residual error R k = I − A S k for k times iteration can be expressed as Then it can be inferred that By comparing (9) and (13), the difference between the convergence speed of the two iteration forms can be found. If the same value of S 0 is selected, S k can be directly obtained in the kth step via the second-order form as (7) while it requires (2 k − 1)-step first-order iteration of (6).

SIMULATION AND RESULTS
The simulation in this paper is performed in CW mode (i.e., ω = 0) and under the assumption of homogenous and approximately infinite medium. The algorithm can also be applied to arbitrary geometries linearized by analytical approximation or finite element method. The measurement geometry for simulations is illustrated in Figure 1. The optical properties of the media are μ s = 10 cm −1 and μ a = 0.03 cm −1 everywhere for both the excitation and emission wavelengths. The original fluorescent yield η is 0.05 cm −1 in the presence of the fluorescent probes. All the simulations were done in Matlab environment (version 7.0.1) on a 2.79 GHz Intel Pentium IV personal computer. The simulated measurement vector U is computergeneralized by the product of coefficient matrix A and the original distribution of X.
In the offline preiteration step, the approximation of A + is obtained by the iteration of (7) with proper iterative number K. However, in the simulation the iterative method is found to have the semiconvergence property. This is probably due to the accumulated round-off error in the computation. So the optimal iteration number should be determined according to experience or prior information about the system. In our simulation, a pretest with a known distribution of fluorescence yield X is performed to choose K for the particular imaging system. And the mean squared error (MSE) between the original X and the reconstructed X is used as a criterion of the reconstructed quality. We investigated how the MSE changed against iteration times for several imaging systems with different sizes (M measurements and N voxels). Figure 2 shows that there is a relative flat segment where MSE changes very slowly before the iterative number begins to rise significantly. So the proper iteration number can be chosen in this iteration number range. Figure 2 is obtained in a noise-free environment. However, it is also found that when noise exists, the MSE rises earlier than in a noise-free system. For different levers of noise, the iteration numbers where MSE rises are different.
With the iterative result S k and the simulated detection U, the distribution of fluorescent yield can be well reconstructed simply by X k = S k U. In our simulation, X k is then modified by including a nonlinear function f to constrain the reconstructed values to [0, 0.05], that is In the simulation, occasions of single-probe as well as multiprobe are reconstructed for several different dataset sizes. It is proved that the distribution of the fluorescent yield η can be well estimated by the fast algorithm ( Figure 3). It can be seen that the algorithm works well when the measurement number is much less then the reconstructed number.
For different imaging subjects, the weighted matrix A may need to be updated, so it would be desirable to know how the inversion time of the preiteration changes with different-sized datasets. According to the results of convergence of the iteration in Figure 2, 60-time iterations are chosen for all of the following datasets in order to compare the reconstructed timescales. The results are shown in Table 1. It can be inferred that for the same number of measurements M, the computing time is approximately proportional to the number of reconstructed voxels N. However, if N remains constant, when M rises to l · M, the computing time will increase to nearly l 2 times of the original.

DISCUSSION AND CONCLUSIONS
With the preiteration method, we have demonstrated reconstruction of fluorescence concentration by using simulation data based on the analytical model with first-order Born approximation. Although in this paper, the fast algorithm is simply demonstrated with the analytical solution for specific medium geometry, it could also be applied to arbitrary geometries, where the model in (1) is discretized by numerical methods [7,10] or the Kirchhoff approximation [5].
In the simulation, a pretest should be done to determine the proper iteration number. A relationship between the convergence property and the dataset size is also obtained through the investigations and it can be found from Figure 2 that in noise-free environment, the number of iterations when the MSE begins to rise mainly depends on the ratio of the voxels number N and the number of measurements M, but not on the absolute value of them. This result will be helpful for the determination of the proper iteration number. For example, the convergence property of large dataset can be predicted from a smaller one with the same N/M. For a system with fixed measurement size, the larger the reconstructed mesh number is, the later will the MSE curve begin to rise.
The computation burden of the second-order iteration is further investigated in our work. It can be inferred from Xiaolei Song et al.   Table 1.This feature should be very suitable for imaging systems where the number of voxels is always much larger than measurement number such as FODT. The results of both convergence property and computation burden indicate that the algorithm is very suitable for imaging systems in which the boundary measurement number is much less than the inside reconstructed voxels. In addition, the reconstructed images in Figure 3 showed that the algorithm works well for these kinds of system. The most promising feature of the algorithm is the rapid reconstruction speed. It significantly accelerates the reconstruction process in the following two aspects. First, when the weighted matrix stays constant or the variance can be ignored, by allowing the time-consuming iteration to be performed offline, it provides great computational facility, which is just a unique matrix vector multiplication. Second, in the preiteration step, it is the second-order iteration form of (7) that exponentially improves the speed of the iterative process, which makes the algorithm feasible in practice and can be finally applied to FODT with datasets of large size. For example, to reconstruct the same quality images with 60 iterations of (7) (the reconstructed images are shown in Figure 3 and the computing time is shown in Table 1), it will cost about 2 60 iterative steps using Landweber method or the firstorder iterative form, requiring days for the reconstruction. So the first-order form is not practical for FODT of large-sized datasets even in the preiteration step. Therefore, the results demonstrate that the time efficiency of both the preiteration process and the online reconstruction is the most important advantage of the algorithm. It will be helpful to promote the development of real-time image reconstruction systems and dynamic monitoring of molecular activity.