A reconstruction approach in wavelet domain for fluorescent molecular tomography via rotated sources illumination

Fluorescent molecular tomography (FMT) aims at reconstructing the spatial map of optical and fluorescence parameters from fluence measurements. Basically, solving large-scale matrix equations is computationally expensive for image reconstruction of FMT. Despite the reconstruction quality can be improved with more sources, it may result in higher computational costs for reconstruction. This article presents a novel method in the wavelet domain with rotated sources illumination. We use the finite element method for the computation of the forward model. The global inverse problem is solved based on wavelet in conjunction with principal component analysis. The iterative reconstruction is implemented with sources rotated in a certain angle. The original excitation light sources are used to reconstruct the image in the first iteration. Then, upon the sources are rotated by a certain angle, they are employed for the next iteration of reconstruction. Simulation results demonstrate that our method can considerably reduce the time taken for the computation of inverse problem in FMT. Furthermore, the approach proposed is also shown to largely outperform the traditional method in terms of the precision of inverse solutions. Our method has the capability to locate the inclusions. The proposed method can significantly speed up the reconstruction process with the high reconstruction quality.

absorption of light at the fluorophore's excitation wavelength. Upon radiative relaxation, fluorescent light is emitted and the fluorophore returns to its ground state with some characteristic time constant. The fluorescent photons are measured by the detectors widely spaced over the surface of the object. From these data, one can detect and map the accumulation of indocyanine green in tissue. Compared to other tomography methods, FMT offers several distinct advantages in terms of sensitivity to functional changes, safety, and cost [6]. For model-based iterative image reconstruction, the light propagation model is utilized as a predictor of measurements. Typically, the model is described by coupled partial differential equations [7]. Besides the forward model, the inversion technique is also needed for image reconstruction [8]. These techniques take into account the diffuse nature of photon propagation to achieve the spatial distribution of fluorochromes in tissues.
Considering the fact that the fluorophore is excited by the excitation light from source, the source may be an important factor for yielding the reconstructed results. Intuitionally, more sources can result in improved reconstruction results. But on the other hand, it may lead to the matrix system with larger scale and hence higher computational costs for reconstruction [9]. A model-order reduction method was proposed in [10] to reduce the computational complexity in the system matrix calculation. However, the transformation matrix needs to be constructed with the basis vectors, which possesses relatively high computational requirements. In [11], an efficient algorithm was proposed to locate and characterize the object, where the B-spline model and appropriate parameterization were utilized to reduce the number of unknowns. However, this method addresses the problem with only one object.
To accelerate inverse problem of FMT, some compression approaches have been proposed [e.g., wavelet transform, principal component analysis (PCA), etc.]. The most important feature of the wavelet transform lies in the fact that most information of the signal is contained in a small number of entries with other entries being very small and therefore can be neglected. PCA is one of the most widely used feature extraction methods, which aims to obtain the most compact representations of the high dimensional data. Some related research has been conducted in the inverse reconstruction. Ducros et al. applied compression techniques to the measurements acquired with structured illuminations [12]. This method is based on the exploitation of the wavelet transform of the measurements acquired after wavelet patterned illuminations. Correia et al. introduced a method with wavelet-based data and solution compression to improve the efficiency of image reconstruction for fluorescence diffuse optical tomography [13]. This approach preserves the resolution of the forward operator and compresses its representation. In [14], Zhang et al. proposed to use PCA to reduce the dimension of the sub weight matrix, and thus to accelerate the reconstruction process of dynamic FMT. Cao et al. solved the inverse problem based on reducing the dimension of the weight matrix with PCA [15]. Furthermore, some other fast reconstruction techniques have been investigated, including sparsity regularization based on the iterated shrinkage method [16], acceleration strategy using graphics processing unit [17,18], and sparsity adaptive subspace pursuit method [5]. In addition, a reconstruction method using permissible region extraction strategy was proposed in [19].
Considering the compression characteristic of wavelet transform and PCA, to further speed up the reconstruction process of FMT as well as improve the precision of the inverse solutions, a new method using the wavelet-based PCA is proposed in this paper. In our method, the original excitation light sources and those rotated in a certain angle are used for iteration of image reconstruction in turn. Simulation results demonstrate that the proposed method can significantly speed up the reconstruction process and achieve high accuracy of inverse solutions.

Diffusion model
As it has been stated earlier, the forward model is used to predict the observable states at the measurement locations from knowledge of the excitation light source and spatial distribution of optical and fluorescent properties. The propagation of photons through a highly scattering medium with low absorption can be well described by the diffusion equation [20]. We employ the widely-used diffusion equation as a forward model that is appropriate for a variety of optical tomography schemes of tissues. Herein, the excitation field � x (r, ω) and the emission field � m (r, ω) are modelled with a pair of coupled diffusion equations as follows where the first equation depicts the transport of the excitation photons and the second one describes the excitation and transport of the fluorescent photons; ∇ is the grad operator, S x (r, ω) is the source term for the excitation light; D x,m (r) and k x,m (r, ω) denote the diffusion and decay coefficients at the excitation and emission wavelengths, respectively; α is the emission source coefficient. They are defined by: where µ ax,mi (r) represent the absorption coefficients due to non-fluorescing chromophore; µ ax,mf (r) represent the absorption coefficients due to fluorophore; µ ′ sx,m (r) denote the isotropic scattering coefficients; fluorescence parameters η and τ (r) denote the fluorescence quantum efficiency and fluorescence lifetime, respectively; c is the speed of light in the media; i is the imaginary unit; ω stands for the angular modulation frequency of the input signal.
Here, we make use of the popular Robin boundary conditions for a bounded domain , which take the form as where n is the outer normal to the boundary, and A x,m (r) is a parameter modelling internal reflection at the boundary.

Finite element approximation of the forward model
Like most others working in FMT, we are currently using the finite element method (FEM) for the computation of the forward model. FEM is versatile especially in regard to complex geometries and for modelling boundary effects [21]. In principle, FEM can be applied to any partial differential equations model of the transport process. In the FEM framework, the computational domain is discretized to a mesh with P elements and N vertex nodes [22]. The solution x,m is approximated by the piecewise function � x,m = N i xi,mi ϕ i , with locally supported basis functions ϕ i (i = 1, 2, . . . , N).
To yield the weak solutions of the forward equations, we rewrite Eqs. (1) and (2) Combining Eqs. (12) and (15), the forward equations within the FEM scheme become where

Inverse problem
The inverse problem of FMT consists in estimating the optical parameters and fluorescent properties of the tissue by using the measured data as described earlier [23]. To generally pose the inverse problem, we first define the forward mapping as F. Therefore, the inverse problem reads where y denotes boundary measurement, and x denotes optical or fluorescent properties. The above non-linear problem can be linearized. To proceed, we can expand about x 0 in a Taylor series. Neglecting the higher order terms, we thus arrive at the linear problem as where J is the Jacobian of the forward mapping.
Due to the fact that the inverse reconstruction problem is ill-posed and underdetermined, we introduce the Moore-Penrose inversion in conjunction with Tikhonov regularization, leading to the following formula: where I represents the identity matrix, and ξ acts as a regularization parameter.
Equation (23) can be written in a succinct matrix form by where we define K = J T J + ξ I and b = J T y.

Image reconstruction with the wavelet-based principal component analysis
We solve the inverse problem in the wavelet domain. To this aim, we take the wavelet transform on both sides of Eq. (24) where However, the level-by-level implementation scheme in the conventional wavelet-based reconstruction method [24] not only is computationally expensive but also causes information lost in the system matrix of the reconstruction problem [25], which inevitably deteriorates the final reconstruction quality. In order to circumvent that problem, we propose to solve the global inverse problem as Eq. (24) based on wavelet in conjunction with the PCA instead of the levelby-level wavelet transform scheme. To this aim, we briefly present the PCA principles. It is well known that PCA performs a dimensionality reduction by searching for a projection matrix with a small number of eigenvectors with respect to the largest eigenvalues. Assume that L is the covariance matrix of the given matrix K, that is, where is a diagonal matrix consisting of the eigenvalues of L, and is the matrix of eigenvectors of L.
Thus the principal components of the matrix K can be achieved by Multiplying (24) from the left with , one has Keeping the first q largest principal components, we can obtain a new matrix equation with reduced scale, namely Therefore, the global matrix system as Eq. (24) can be approximately solved with the reduced-scale matrix system according to PCA.
The inverse reconstruction with the wavelet-based PCA is summarized in "Algorithm 1".
1. Take wavelet transform with respect to K and b in Eq. (24) to achieve the approximation components K 1 and b 1 ; 2. Solve K 1 x 1 =b 1 with PCA; 3. Prolongate x 1 by padding zeros to achieve an initial guess for x at the original resolution, i.e.,

Iteration based on the strategy of excitation light sources rotation
The tomographic imaging involves placing sources and detectors over the available surface of the tissue. Basically, the excitation light sources are arranged at the fixed positions during the process of image reconstruction. By means of increased number of sources, the image quality can be improved. However, such a strategy may result in the matrix system with larger scale and hence higher computational complexity. Although we can reduce the number of sources to safe the computation time, the information for image reconstruction will decrease, which may lead to the poor quality of reconstruction. As a result, there exists a contradiction between the reconstruction accuracy and the computational requirements. In order to address this tradeoff, we propose a new strategy for iterative calculation. In such a strategy, the original excitation light sources are used to reconstruct the image in the first iteration. Then, upon the sources are rotated by a certain angle, they are employed for the second iteration of reconstruction. This means that the whole iterative reconstruction is performed using the sources with different rotation angles in turn. This process is repeated until some stopping criteria are satisfied. This strategy is motivated by the fact that the excitation light sources from different angles can provide more information than those from some fixed angle during the iteration process, and thus the quality of reconstructed results can be improved. In our method, the number of excitation light sources is not increased, and thus it will not lead to higher computational cost. Moreover, the iterative results from the sources with one angle can provide a good initial guess for the next iteration from the sources with other angle. In this way, the precision of solutions can be improved with rotation of the lights. However, if the rotation angle of sources is too small, it may provide quite limit information for reconstruction. On the contrary, large rotation angle may lead to the superposition between the original and rotated sources. In such a case, it is unable to provide additional information for iterative reconstruction. Suppose the sources are distributed around the circumference of the tissue with equal angle between each source. To overcome those difficulties, in our work, the rotation angle is set as a half of the angle between each source. This strategy can be schematically illustrated as in Fig. 1.
For derivation of the algorithm, we minimize the residual error between the predicted data and measured data to acquire the solution to the reconstruction problem by where M(x) is the objective function, y is the measured data, and F (x) is the predicted data with regard to a forward model. Let us suppose that β is a half of the angle between each source, and thus the resulting reconstruction algorithm is summarized in "Algorithm 2".

Simulation results and discussion
In this section we performed simulation study using different phantoms to test the performance of our algorithm and the obtained results. The forward model as Eqs. (1) and (2) is used to simulate the measured data. In order to better simulate the realistic conditions, we add Gaussian noise with a signal-to-noise ratio of 10 dB to the calculated data. Actually, large regularization parameter may lead to low contrast and resolution of the image, while small parameter can result in increased contrast and resolution. However, small parameter also increases the high frequency noise in the image [26]. The regularization parameter ξ is set to 0.001 in the simulations for better results after a lot of simulations. The termination criterion δ is set to 0.02.
In the first example, verification of the performance of the proposed method is investigated using the test phantom containing one inclusion as indicated in Fig. 2. Four excitation light sources are uniformly distributed around the simulated phantom. The measurements are sampled by thirty detectors uniformly placed on the boundary of the phantom.
To reduce the computational requirements without significant reduction of image resolution, we compute the reconstructions based on the mesh that is adaptively refined with respect to the a priori image as portrayed in Fig. 3. Figure 4 displays the mesh containing 122 nodes and 212 triangular elements. The details for optical and fluorescent parameters in different areas of the test phantom are provided in Table 1. In order to compare the reconstructed object with the true one, we define an image quality metric by introducing the mean square error (MSE), given as   where the superscript rec denotes the values obtained using reconstruction algorithms, and act denotes the actual distribution of the optical or fluorescent parameters which is used to generate the synthetic image data set.
The reconstructed images of µ axf for one-inclusion phantom with two sources and four sources are depicted in Fig. 5a, b, respectively. Both of them are obtained without using the wavelet-based PCA. The results presented in Fig. 5 can be explained by considering that reconstruction with the increasing number of sources can enhance the quality of image, whereas the time requirements for reconstruction may increase.
In Fig. 6 we show the resulting reconstructions using the different algorithms. Figure 6a displays the reconstructed result using the proposed method with four sources. Figure 6b, c depict the traditional reconstructed result with four sources and that with eight sources, respectively. We see that the method proposed is capable of yielding the reconstructed target with improved contrast and contour comparatively to the traditional method.
We demonstrate the benefits of the proposed method by comparing the performance of our method to the traditional method. For quantitative validation, the performance of reconstructions in terms of the computation time and MSE is tabulated in Table 2. We remark, that the computation time for the proposed algorithm is much faster than the traditional method, which demonstrates that our method is time efficient. Although the increased number of sources can improve the quality of reconstruction, it will slow down the speed of reconstruction. In addition, the MSE of the proposed method is Table 1 Optical parameters used for one-inclusion phantom smaller than that of the other compared method. Therefore, the above results suggest that the algorithm proposed can substantially speed up the reconstruction process and possess high accuracy. The phantom for the second test case is shown in Fig. 7. It consists of two inclusions of different shapes. As before, the phantom is illuminated by four equally spaced sources located on its boundary. The detector readings are obtained from 30 different points from the boundary of the circular domain. The distance between the successive detector positions is the same through the boundary. Figure 8 displays the a priori image as a guidance for generation of the adaptively refined mesh. The resulting mesh with 148 nodes and 264 triangular elements is depicted in Fig. 9. Table 3 lists the values of optical and fluorescent parameters of the simulated phantom.   In Fig. 10 we show the reconstructed images of µ axf for two-inclusion phantom with 2 sources (see Fig. 10a) and that with four sources (see Fig. 10b). We also notice that one can obtain better reconstructed results with increasing sources. Nevertheless, reconstruction with more sources may lead to a heavy computation burden.
The reconstruction from our method with four sources is shown in Fig. 11a and those from the traditional method are depicted in Fig. 11b, c. Particularly, the reconstructed images are obtained with four sources (see Fig. 11b) and eight sources (see Fig. 11c).
We find that the contrast of image can be enhanced with our method. The reconstruction with more sources can improve the reconstruction accuracy while result in a heavy   Table 3 Optical parameters used for two-inclusion phantom computation burden. More importantly, we also note that the proposed method can improve the quality of reconstruction with more accurate shape and position of both targets. We provide the quantitative comparisons of different reconstructions presented in Table 4. As can be clearly seen, improvement in quality of reconstruction can be achieved by proposed algorithm. Additionally, it is evident from Table 4 that our method requires less reconstruction time as compared with the traditional method. Therefore, the main conclusion we can draw from these simulation studies is that the approach proposed has comparable computational efficiency to the traditional method and high capability to achieve accurate reconstruction.   To illustrate the superiorities of the proposed algorithm, we show the reconstructed results from the different algorithms (see Fig. 12). Figure 12a-c display the reconstructed images with the proposed approach, wavelet method, and PCA method, respectively. Table 5 summarizes the quantitative performance of reconstruction. From Table 5, it can be clearly seen that the proposed algorithm has better performance on accuracy and speed of reconstruction than algorithms only using wavelet method or PCA method.
To validate the proposed approach in the 3D case, the methods previously defined for triangular elements are extended to tetrahedral elements. The integration of products of shape functions over the volume of the elements, and surface integrals over a side of the element is performed by numerical integration rules. Herein, a cylindrical phantom as illustrated in Fig. 13 is utilized for 3D simulations. A small cylindrical inclusion is suspended in this phantom. The dashed curves represent the planes of measurement. Six sources and sixteen measurements are employed for each plane. The data are collected in all three measurement planes. The mesh for 3D reconstruction containing 3208 tetrahedral elements as well as 858 nodes is shown in Fig. 14. Figures 15 and 16 display the 3D reconstructed images based on the proposed approach and the traditional method, respectively. These are 2D cross sections through the reconstructed 3D images. The quantitative performance of the above two methods is given in Table 6 to further evaluate the reconstruction quality. As one can see from Table 6, our proposed algorithm can also significantly speed up the reconstruction process and improve the quality of reconstruction in the 3D case.
Finally, we test the reconstruction algorithms with the Monte Carlo method. As most commonly used stochastic technique, Monte Carlo method is regarded as gold standard for modelling the light propagation and has a long pedigree in transport theory. We utilize the Monte Carlo method to generate the measurement data, which is employed to reconstruct the image of FMT. Figure 17 shows the model for reconstruction and Fig. 18     shows the corresponding reconstructed results with four sources. The reconstructed results are obtained from the proposed algorithm (see Fig. 18a) and the conventional method (see Fig. 18b). The quantitative performance is listed in Table 7, from which we can also see that both the speed and precision of the reconstruction can be improved with the proposed algorithm.