Abstract
Image reconstruction of fluorescent molecular tomography (FMT) often involves repeatedly solving large-dimensional matrix equations, which are computationally expensive, especially for the case where there are large deviations in the optical properties between the target and the reference medium. In this paper, a wavelet-based multiresolution reconstruction approach is proposed for the FMT reconstruction in combination with a parallel forward computing strategy, in which both the forward and the inverse problems of FMT are solved in the wavelet domain. Simulation results demonstrate that the proposed approach can significantly speed up the reconstruction process and improve the image quality of FMT.
1. Introduction
The use of near-infrared (NIR) light in biomedical research has made significant progress over the past few years [1]. It has been shown that light with wavelengths in the near-infrared range can propagate through tissue for distances on the order of multiple centimeters, because of low tissue absorption in the “near-infrared window.” This finding has encouraged the development of fluorescence techniques to visualize specific biochemical events inside living subjects (in vivo molecular imaging) [2]. Fluorescence techniques have played a critical role in the description of biological processes at the molecular and cellular levels [3]. One particular example is fluorescent molecular tomography (FMT), which is an emerging tool for molecularly based medical imaging [1]. In this imaging modality, a fluorescent biochemical marker used as contrast agent is injected into the biological system and consequently accumulates in diseased tissue as a result of the increased vascular density or by means of selective targeting [4]. During the imaging process, light at the fluorophore’s excitation wavelength is launched into the tissue, and then it is absorbed by fluorophore that presents in the tissue, and the fluorophore is elevated to an excited state and remains there for some period of time (the fluorescence lifetime). Some proportion of the excited molecules will ultimately release their excess energy by emitting a photon as they drop back to the ground state. This creates fluorescence which can be separated from the excitation light via interference filters [5]. Volume images of the fluorescent yield and lifetime parameters are reconstructed from several optical measurements on the surface of the tissue [4].
Reconstruction of tomographic data from diffusing sources involves the generation of a forward model that predicts the photon distribution striking the detectors for a given source location and medium [6]. One challenging problem in the reconstruction process is that the computational complexity is very high due to an extremely large dimension of the matrix, which is not only in the inverse problem but also in the forward problem. Multiresolution approach is an effective way to speed up the process of solving the above problem. It is well known that the most important feature of the wavelet transforms 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. In [7], an efficient pyramidal algorithm was proposed for the multiresolution representation of the signal with wavelets using orthogonal basis functions and quadrature mirror filters to compute it. Unser et al. extended these ideas to the case of nonorthogonal basis functions using splines [8]. Because the forward problem should be repeatedly solved during the process of solving the inverse problem especially for the case where there are large deviations in the optical properties between the target and the reference medium, the speed and accuracy of the forward computation are of critical importance determining the performance of the reconstruction algorithm. In order to speed up the forward computation process, we propose to generalize the strategy in [9] for solving the forward problem of FMT in the wavelet domain in combination with a parallel computing strategy [10]. Simulation results demonstrated that the proposed algorithm can significantly improve the efficiency of reconstruction for FMT. The main contribution of this paper is the extension of the multiresolution reconstruction approach originally developed for the diffuse optical tomographic reconstruction to the case of fluorescent molecular tomographic reconstruction suitable for the case where there are large deviations in the optical properties between the target and the reference medium. The forward problem of FMT is solved in wavelet domain in combination with a parallel computing strategy originally developed by our group, which decouple the two originally coupled differential equations corresponding to the excitation and the emission light, making the forward problem suitable for parallel implementation [10].
2. Forward Problem
2.1. Governing Equations
The generation and propagation of the excitation and fluorescence (emission) light in random highly scattering media can be described by two coupled diffusion equations which are the P1 approximation to the radiative transport equation (RTE). In the frequency domain, the diffusion equations become elliptic and can be expressed as subjecting to the Robin boundary conditions on the boundary of the tissue where (1) describes the propagation of the excitation light, and (2) models the generation and propagation of fluorescently emitted light. The subscripts and denote the excitation and emission light wavelengths, respectively. is the gradient operator, is the excitation light source, and is the photon fluence. is a vector normal to the boundary , and are the Robin boundary coefficients which are governed by the reflection coefficients (,). The values of and are 1/2 (no reflection). In addition, the diffusion coefficients , decay coefficients , and emission source coefficients are, respectively, defined as where and are the absorption coefficients due to nonfluorescing chromophore; and are the absorption coefficients due to fluorophore; and are the isotropic scattering coefficients; is the fluorescence quantum efficiency; is the fluorescence lifetime; is the speed of light in the media and .
2.2. Finite Element Formulation
The solutions to (1) and (2) can be obtained using the finite element method (FEM) which is a completely general numerical technique applied to any geometry [11]. The FEM is one of the most popular methods for numerically solving partial differential equations (PDEs) because of its applicability to a range of problems and the existing large body of mathematical theory [12]. In the FEM framework, the domain is divided into elements, joined at N vertex nodes. The solution is approximated by the piecewise function , with being basis functions [13].
Suppose [14], , we have
Let , where stands for both and . In order to obtain the weak solutions of (1) and (2) under the boundary conditions of (3), (1) and (2) are written as the following sesquilinear form: where
where and are, respectively, the bounded domain and its boundary. Equation (6) can be further rewritten in a more compact matrix form as where The elements of finite element matrix can be obtained from the following formula:
2.3. Forward Computation
The accuracy and speed of solving the forward problem as discussed in Section 1 are of critical importance determining the performance of the reconstruction algorithm. A multiresolution iterative perturbation reconstruction method for the optical tomographic image reconstruction based on the wavelet transform is presented in [9]. By computing the Jacobian matrix, which is a measure of the rate of changes in measurement with respect to the optical parameters, at a reference medium whose optical properties are similar to those of the target medium, the reconstruction problem is reduced to a system of linear equations. As a result, there is no need for repeated solving of the forward differential equations. However, this method would not be valid for the case where there are large deviations in the optical properties between the target and the reference medium. In such a case, the forward problem should be repeatedly solved during the process of solving the inverse problem. In order to speed up the forward computing process so as to speed up the whole process of reconstruction, we propose to generalize the strategy proposed in [9] originally developed for the inverse reconstruction of the diffuse optical tomography to the case of the forward problem of FMT and solve it in the wavelet domain. Furthermore, in order to decouple the forward problem of FMT, a parallel strategy previously developed by our group [10] will be used in combination with the aforementioned strategy for solving the forward problem. Our innovations are especially suitable for the case where there are large deviations in the optical properties between the target and the reference medium.
2.3.1. Brief Introduction of Wavelets
For the convenience of the following discussion, a brief introduction of the theory of wavelet transform is presented here. The wavelet transform is a tool that cups up data of functions or operators into different frequency components and then studies each component with a resolution matched to its scale. By a proper design of the basis, the wavelet can project the signal onto a chain of embedded approximations and details at various levels of resolutions, and, as a result, the wavelet transform is usually referred to as the multiresolution analysis. For example, the two-level wavelet-based multiresolution representation of one dimensional discrete signal with components can be described as where is the wavelet transform of the original signal . It can be seen from this equation that the original signal can be decomposed into two parts of the detail component and the approximation component .
Similarly, the two-level wavelet-based multiresolution representation of a 2D image sized can be expressed with the following formula: Four elements in the matrix of the right-hand side of (13) are, respectively, the approximation image and three detail images , , and .
2.3.2. Multiresolution Computing of the Forward Problem for the Excitation Light in Wavelet Domain
In order to exploit the multiresolution property of the wavelet and reduce the forward computational time, the forward problem is first represented in the wavelet domain. For such a purpose, multiplying both sides of (8) from the left by and assuming the orthonormality of , we have where , , , and are, respectively, the wavelet transform matrix of and Φ.
It is well known that the most important feature of the wavelet transforms 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. As a result, the dimension of the forward problem can be reduced level by level by using only the approximation components of the wavelet coefficients to describe the forward problem, that is, where denotes the index of the scale, , , , and and are, respectively, the wavelet transform matrix of and at the th scale with and being, respectively, the approximation components of the corresponding signal at the th scale, and is the LL components of the corresponding wavelet transformed stiffness matrix at the th scale.
Using the above multiresolution representation, the forward problem can be solved in a fine-to-coarse-to-fine procedure which can be summarized as in Algorithm 1.
| ||||||||||||||||||||||||||||
Owing to the fact that some important features are contained in the coarse resolution solution, as a result, it will be very helpful for speeding up the iterative process when solving the forward problem at a higher level resolution with the solution obtained at a coarser resolution as an initial guess. Therefore, we can expect to expedite the process of solving the forward problem by using Algorithm 1 with a fine-to-coarse-to-fine strategy.
2.3.3. Parallel Implementation of the Forward Problem for the Emission Light
After the discussion of the wavelet-based algorithm for the forward problem corresponding to the excitation light, the next task for us will be that of solving the forward problem for the emission photons. For the case where there are large deviations between the referenced and target medium, the forward equations must be solved repeatedly during the process of reconstruction following a model-based iterative image reconstruction scheme. Therefore, a rapid and accurate computational implementation of the forward problem is of critical importance for fluorescent molecular tomographic image reconstruction. From (8) and (9), we can see that the two forward models corresponding to the excitation and emission light at different wavelength coupled together because the solution to (8) is contained in the source term of (9). Traditionally, the forward problem of (8) and (9) are solved in a sequential manner, that is, (8) is first solved whose solution is then substituted to (9), which yields the photon fluence at the emission wavelength. That scheme will affect the computational speed of the forward problem, even the inverse problem. To tackle such a problem, an approximate computing strategy for decoupling these two forward equations was proposed in [15] and was used for the FMT reconstruction in [16]. However, this strategy is not valid for the case where there is a large stokes shift [15]. For a rapid implementation of the forward problem, we have proposed an accurate parallel implementation scheme in [10] where the following equation is solved instead of (9): In (16), is an identity matrix. Since is symmetric and positive definite [10], we can always obtain an inverse matrix for from (16). The matrix can be obtained with the numerical method which can be speeded up when the matrix is symmetric and positive definite [17]. Combining (8) and (16) leads to a system of equations in discretized domain for the forward problem of FMT. Because (8) and (16) are independent, they can be solved in a parallel manner. Obviously, the photon fluence of at the emission wavelength can be recovered by simple matrix multiplication with contained in obtained from Algorithm 1, that is,
In summary, the whole forward computation process in our proposed algorithm can be realized with Algorithm 1 in the wavelet domain for the excitation light and in a parallel manner for the emission light. It has been proved in our simulations that our proposed forward computing algorithm can significantly reduce the computational requirements.
2.3.4. Computational Complexity Analysis of the Parallel Computing Strategy
The most important aspect of the parallel computing strategy is decoupling of the two coupled equations. In order to illustrate the improvement of the parallel computing strategy in computational complexity as compared with the sequential one quantitatively, the computing efficiency is analyzed as follows.
Because the maximum computational complexity of solving linear equations defined by a matrix sized with CGD method is ) (the maximum iteration number for such an optimization problem is [18], and the complexity of each iteration is ) [19], the complexity of solving the coupled forward problem of FMT, that is, (8) and (9), in a sequential manner will be . On the other hand, if the forward problem is solved according to the parallel strategy as discussed previously, the computational complexity will be because the operations of matrix inversion in (16) and solving the linear equation in (8) can be implemented independently with two processors simultaneously whose computational complexities are , while the computational complexity of multiplying the matrix by a vector in (17) is [20]. Thus, the speed of forward computing can be improved in such a parallel manner. The above analysis is valid for both two-dimensional (2D) and three-dimensional (3D) cases, because the only difference between these two cases is that the size of the matrix for 3D case is much larger than that of 2D case.
Particularly, if we are interested only in the reconstruction of the absorption coefficient due to the fluorophore, the matrix in (16) needs to be calculated only once during the whole reconstruction process and hence the computational requirements are extremely reduced.
3. Image Reconstruction of FMT
3.1. Inverse Problem
The forward and inverse problem of FMT can be, respectively, formulated as where is the detector readings, is the forward operator, and x is the optical or fluorescent properties of the tissue.
Generally, is a nonlinear function of . In order to simplify the reconstruction process, we expand the function in the vicinity of in a Taylor series [21]:
where and are the first- and second-order Frechet derivatives of and are usually referred to as the Jacobian matrix and Hessian matrix, respectively, if represented in matrix form. Keeping up to the first-order terms in (19) and introducing the Tikhonov regularization term for tackling the ill-posedness of the inverse problem, the linearized formulation for the reconstruction problem can be described by where is an identity matrix, is a regularization parameter, and is the Jacobian matrix describing the influence of a voxel on a detector reading [22]. The Jacobian matrix is obtained using the perturbation method which can be described as
where is the perturbation in the optical or fluorescent properties, and stands for the corresponding residual data between the two predicted data.
By introducing two quantities of and which are, respectively, the perturbation in the optical or fluorescent properties and residual data between the measurements and the predicted data, (20) can be rewritten in a more compact matrix form as [23] with and
Using (22), we can obtain the reconstructed image simply by finding a solution of to it. In our case, both the matrix and the vector and hence the Jacobian matrix J and the residual data are functions of considering the fact that there are large deviations between the target and reference medium, which is different from the case in [9] where and are irrelative to and remain unchanged during the iteration process. Therefore, both and should be repeatedly calculated during the reconstruction process if iterative method is used to find a solution to (22), which means that the forward problem should also be repeatedly solved in the reconstruction process. As a result, the reconstruction efficiency can be significantly improved if we can expedite the repeated forward computation. As mentioned before, we can expect to expedite the process of the forward computing with an algorithm in the wavelet domain discussed in Section 2. In order to further speed up the reconstruction algorithm, we propose to adopt the multiresolution reconstruction scheme in the inversion process. For such a purpose, we perform the wavelet transform on both sides of (22) and have
where , , , and are, respectively, the wavelet transform matrix of and , and is an orthonormal matrix. As for the case of the forward problem, the wavelet transform can be successively performed level by level with respect to the approximation components of both sides of (23) and obtain a multiresolution representation of the reconstruction problem. The whole reconstruction algorithm can be summarized as in Algorithm 2.
| ||||||||||||||||||||||||||||||||||||||||
From Algorithm 2, we can see that there are actually two layers of iterations in it: one is the inner iteration where the Jacobian matrix is not updated which is similar to that proposed in [9], and the other is the outer iteration where both the Jacobian matrix and the residual vector are recomputed at the new values of the optical parameters obtained in the former iteration. Owing to the fact that both the Jacobian matrix and the residual vector are updated during the outer iterations, our algorithm can free the constraints of small deviations of the properties between the target and the reference medium.
3.2. Data Correction
Actually the fluorescent measurements are used as the input to reconstruct the image for FMT according to Section 2. Usually the fluorescence may exist not only in the target but also in the background [24]. When the fluorescence image is reconstructed, it may contain the target fluorescence as well as the background fluorescence. Therefore, if the detector readings are directly used for image reconstruction, the performance of the reconstruction result will drop [25]. In order to improve the reconstruction quality in the presence of the background fluorescence, the data of reconstruction need to be corrected.
In the presence of the background fluorescence, the fluorescence concentration can be formulated as follows [25]: where is the total fluorescence concentration, and and denote the target and background fluorescence concentration, respectively.
Furthermore, the fluorescence concentration can be described as follows [25]: where is the fluorescence quantum efficiency, and is the absorption coefficients.
According to (24) and (25), the absorption coefficients of the fluorescence can be formulated as the sum of the absorption coefficients of the target and background fluorescence, that is,
where and are the absorption coefficients of the target and background fluorescence, respectively.
In order to improve the reconstruction quality in the presence of the background fluorescence, the reconstruction results can be corrected as follows [25]:
where is the domain of interest for the reconstruction with an area of . From (27), it can be seen that can be obtained by taking the average of the reconstructed absorption coefficients especially when the variations of background fluorescence concentration are small. Additionally, as no negative fluorescence exits, any negative value of should be set to zero.
4. Simulation Results
4.1. Two-Dimensional Reconstruction
The algorithm proposed in this paper has been firstly tested in a 2D simulated phantom with two anomalies existing within it as illustrated in Figure 1. Four sources and thirty detectors equally distributed around the circumferences of the phantom are adopted in the simulations. The optical and fluorescent parameters in different regions of simulated phantom are listed in Table 1. The simulated forward data are obtained from (1) and (2), in which the Gaussian noise with a Signal-to-Noise Ratio of 15 dB is added for evaluating the noise robustness of the algorithms. Furthermore, the background fluorescence is also included in the simulated data for evaluating the data correction strategy according to Table 1. Since Daubechies 1 (haar wavelet) has the advantages such as orthogonality and symmetry, Daubechies 1 wavelet as defined in (28) is used in the simulations [26]:
In our current implementation, we will focus on reconstructing the distributions of absorption coefficients due to the fluorophore. The termination criterion in Algorithm 2 is set to 0.02. The regularization parameter is set to 0.005 in the simulations for better results after a lot of simulations [27]. The uniform mesh for reconstruction is shown in Figure 2 with 91 vertex nodes in it. Two quantities are introduced for the quantitative evaluations of different algorithms. The first one is the error function between the simulated data and the predicted data computed at the final reconstructed value, that is,
where is the number of vertex nodes, is the detector readings, is the forward operator, is the reconstructed result of optical or fluorescent properties of tissue, and is -norm. The second one is the normalized root mean squares (NRMSs) error for the reconstructed results defined as where is the number of vertex nodes, and are the original pixel and reconstructed pixel values, respectively, and is the mean value of the original pixel.
The data correction strategy is implemented after the reconstruction for improving the reconstruction quality. Figure 3(a) depicts the reconstructed result of absorption coefficients without data correction, and Figure 3(b) shows the corresponding result after data correction, both of them are based on the proposed wavelet-based multiresolution algorithm. From these two images, it can be seen that the data correction strategy can improve the image quality. Hence, all reconstruction results presented in the latter part of this section are those after data correction.
(a)
(b)
Figures 4(a) and 4(b) show the reconstructed images of using the proposed algorithm and the method in [9], respectively. In this example, the deviation of the optical properties between the reference and the target medium is set to a relatively larger value (here is set to 10 ) for an illustration of the reliabilities of the above two different algorithms under the circumstance with large deviations. These two algorithms are implemented with a same initial guess of 10. Table 2 summarizes the performance of these two algorithms in terms of and NRMS. From this table it can be seen that the proposed algorithm can achieve a more accurate reconstructed result for the case where there are large deviations in the optical properties between the target and the reference medium, that is, our algorithm is more suitable for such a case.
(a)
(b)
Recently, there has been a great amount of interest in developing multimodality imaging techniques for oncologic research and clinical studies with the aim of obtaining complementary information and, thus, improving the detection and characterization of tumors [28]. As a result, it will be helpful to incorporate the prior information obtained from other imaging modalities in the reconstruction process for reducing the computational requirements while achieving a relatively better reconstructed result. In our case, we proposed to use the prior information to generate a nonuniform mesh for the fluorescent image reconstruction according to an adaptively refinement scheme. The basic idea of this scheme is that, for areas in the prior image with large variations of the pixel values, the mesh at this position should be locally refined, and hence the image will be reconstructed with higher resolutions, whereas for regions with small variations, the mesh at this position should be left unchanged, and hence the image will be reconstructed with low resolutions locally at this position correspondingly. Obviously, this idea is plausible because flat regions contain little information and therefore low resolution reconstruction will not lead to serious degradation of the reconstructed results. However, this nonuniform reconstruction will significantly reduce the computational requirements as compared with the uniform fine reconstruction. To simulate such an idea in fluorescent image reconstruction, we use the image shown in Figure 5 with a resolution of pixels as a prior image. Figure 6 shows the adaptively refined mesh with 148 vertex nodes in it which is generated based on the above idea. All of the following reconstructed results of this section are obtained based on this adaptively refined mesh.
Furthermore, in order to demonstrate the advantage of the proposed algorithm as compared with the traditional single resolution method without wavelet transform, Figures 7(a) and 7(b) depict the reconstructed images using the proposed wavelet-based multiresolution reconstruction algorithm and the single resolution method, respectively, both of which are implemented in combination with the parallel forward computing strategy for further speeding up the process of the reconstruction. Table 3 summarizes the performance of the different algorithms in terms of the computation time, NRMS and , from which we can see that our wavelet-based algorithm outperforms the single resolution algorithm in both the reconstruction accuracy and the computational requirements.
(a)
(b)
4.2. Three-Dimensional Reconstruction
To further validate the proposed algorithm for 3D reconstruction, we extend the methods previously defined for triangular elements to tetrahedral elements. Therefore, the shape functions in the local coordinate system is defined as
The integration of products of shape functions over the volume of the elements, and surface integrals over a side of the element, as required for the computation of element stiffness and mass matrices, is performed by a numerical integration rules. Once the element matrices are computed, the FEM model can be solved as in the 2D case without needs for any further alteration.
In the 3D case, a phantom of radius 10 mm and height 40 mm with a uniform background as illustrated schematically in Figure 8 is used for simulations. Within this phantom, a small cylindrical object of radius 2 mm and height 6 mm with is suspended. In Figure 8, the dashed curves represent the planes of measurement, 5 mm apart with coordinates of 15, 20, and 25 mm. Four sources and sixteen measurements are used for each plane in the simulations. The mesh for reconstructing the 3D image as shown in Figure 9 is a cylindrical mesh of radius of 10 mm and height 40 mm. It contains 858 nodes and 3208 tetrahedral elements. The data are collected in all three measurement planes, as shown in Figure 8. Figures 10 and 11 depict the 3D reconstructed images using the proposed algorithm and the single resolution method, respectively. These are 2D cross sections through the reconstructed 3D images. The right-hand side corresponds to the top of the cylinder (), and the left corresponds to the bottom of the cylinder (), with each slice corresponding to a 10 mm increment in the z coordinates.
Table 4 lists the performance of the above two methods for a quantitative comparison in detail. It can be seen that the proposed algorithm can significantly speed up the process of reconstruction and improve the reconstructed image quality. Therefore we can conclude that our proposed algorithm also outperforms the single resolution reconstruction algorithm without wavelet transform for the 3D case.
Furthermore, the proposed algorithm decouples the two coupled equations for the forward problem of FMT, and thus it is quite suitable for parallel computing of the two independent equations with two processors. Table 5 summarizes the computation time and rate of speedup with different number of processors used to validate the superiority of the proposed parallel computing strategy. From Table 5, it can be seen that the parallel computing strategy can speed up the reconstruction process both in the 2D and 3D reconstruction. The rate of speedup for the 3D case is a little higher than that for the 2D case, which indicates that the superiority of parallel computing strategy is more prominent in 3D reconstruction than in 2D case.
5. Conclusion
In summary, a wavelet-based multiresolution reconstruction algorithm is proposed in combination with the parallel forward computation strategy for the purpose of speeding up the reconstruction process with an improved reconstruction accuracy. The most important contribution of this paper is the novel extension of the multiresolution reconstruction approach originally developed for the diffuse optical tomographic reconstruction to the case of fluorescent molecular tomographic reconstruction and for the case where there are large deviations of the optical parameters between the target and the reference medium. Different from the algorithm proposed in [9], the forward problem of FMT is solved in wavelet domain in combination with a parallel computing strategy for speeding up the forward computing process which is especially suitable for the case where there are large deviations in the optical properties between the target and the reference medium, and thus the forward problem should be computed repeatedly. Simulation results demonstrate that the proposed algorithm can significantly reduce the computational complexity and achieve a higher reconstruction quality.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (no. 60871086), Natural Science Foundation of Jiangsu Province China (no. BK2008159), ARC, and PolyU Grants.