Fluorescence molecular tomography (FMT), as a promising imaging modality, can three-dimensionally locate the specific tumor position in small animals. However, it remains challenging for effective and robust reconstruction of fluorescent probe distribution in animals. In this paper, we present a novel method based on sparsity adaptive subspace pursuit (SASP) for FMT reconstruction. Some innovative strategies including subspace projection, the bottom-up sparsity adaptive approach, and backtracking technique are associated with the SASP method, which guarantees the accuracy, efficiency, and robustness for FMT reconstruction. Three numerical experiments based on a mouse-mimicking heterogeneous phantom have been performed to validate the feasibility of the SASP method. The results show that the proposed SASP method can achieve satisfactory source localization with a bias less than 1mm; the efficiency of the method is much faster than mainstream reconstruction methods; and this approach is robust even under quite ill-posed condition. Furthermore, we have applied this method to an in vivo mouse model, and the results demonstrate the feasibility of the practical FMT application with the SASP method.
© 2014 Optical Society of America
With the solid development of fluorescent probes and reporter technologies , the application of fluorescence molecular imaging in biomedical study has become a hot spot over the past few years [2, 3]. It offers an opportunity for noninvasive visualization of biological processes at the cellular and molecular levels. Fluorescence molecular tomography (FMT) is a new promising imaging modality, which can three-dimensionally detect fluorescence biodistribution in vivo, therefore greatly facilitate its applications in small-animal research and pre-clinical diagnostics [4, 5].
However, it is challenging for FMT to reconstruct the fluorescence biodistribution effectively and robustly , because FMT presents a challenging inverse problem which is quite ill-posed or ill-conditioned due to the following three reasons [6, 7]:
- 1. Photons emitted from the fluorescent sources undergo multiple scattering and potential absorption in biological tissues.
- 2. Only the photon distribution on the surface is measurable. Although increase of the measurement data sets can reduce the ill-posedness, the problem may still be ill-conditioned as it is too sensitive to noise caused by charge-coupled-device (CCD) measurement errors and data discretization errors.
- 3. Furthermore, the high sampling measurements and the real animal-shape geometry modeling usually produce a large amount of data, which lead to large computational complexity and big storage capacity.
Therefore, the FMT reconstruction faces various challenges on its efficiency and robustness, and the development of feasible FMT reconstruction approaches plays a unique role in the achievement of practical biomedical applications.
In order to overcome the above challenging problems, some researchers have proposed various regularization methods to make the solution stable and insensitive to noise. Among different regularization methods, the Tikhonov regularization is one of the most popular regularization methods and has been widely applied in resolving FMT problems [8, 9]. It adds an L2-norm constraint of the solution to the original problem, and aims to reduce high-frequency noise in the reconstructed results. The primary benefit of using Tikhonov regularization is that the optimization problem can be simplified and efficiently solved by using standard minimization tools, such as the conjugate gradient method and Newton method. However, it tends to produce an over-smooth solution and loss of some localized features during the reconstruction process by using Tikhonov regularization .
In order to improve the quality of the reconstruction for FMT, more a priori information should be included [6, 11]. Over the past few years, considerable attention has been focused on the sparse regularization being used in the field of compressed sensing (CS) for signal recovery and image processing. According to the CS theory, a sparse or compressive signal can be faithfully recovered from far fewer samples or measurements . Considering that in the practical applications of FMT, the fluorescent sources are usually sparse because the fluorescent probes used in FMT are designed to locate the specific areas of interest such as tumors or cancerous tissues, which are usually small and sparse compared to the entire reconstruction domain . Hence, the problem of FMT can be regarded as a sparse reconstruction problem and the fluorescent source distribution can be recovered by using sparse-type regularization (L0-norm and L1-norm) strategies. Inspired by the ideas behind the CS theory, several algorithms incorporated with L1-norm regularization have been proposed to solve the optical tomography problems in recent years [10, 13–17]. To preserve the sparsity of the fluorescent sources, an iteratively reweighted scheme based approach, which was able to obtain more reasonable and satisfactory results compared with the Tikhonov method was proposed . At the same time, to improve the reconstruction accuracy, an effective FMT reconstruction algorithm based on the iterated shrinkage method with the L1-norm (IS_L1) was proposed , the reconstruction algorithm was able to comparatively acquire accurate results even with quite limited measurement data sets. However, the approximate convergence rate of this algorithm is linear since it is a first-order method . Therefore, it needs a large number of iterations to reach an acceptable solution, especially when the dimension of the FMT inverse problem is quite large. To enhance the reconstruction efficiency, the stagewise orthogonal matching pursuit (StOMP) based method was introduced to conserve computation time due to the greedy pursuit strategy , while it needs to estimate the sparsity factor which indicates the number of unknowns in advance. There are different sparsity factors for different FMT experiments, therefore it is not always able to achieve acceptable results for estimation of the sparsity factor empirically in this method. Although the aforementioned reconstruction methods usually work well in some specific and highly controlled situation, further study is urgently needed to investigate more general cases .
In this paper, a novel method based on the sparsity adaptive subspace pursuit (SASP) has been proposed for FMT reconstruction. This novel method adopts a subspace projection and correlation maximization approach to simplify the FMT problem with sparsity-promoting L1-norm regularization and to treat it as the basis pursuit problem . The proposed method performs reconstruction by employing a search strategy in which a number of K (i.e., the sparsity factor which indicates the number of unknowns) vectors with the highest correlation are selected from the candidate set. Then the search strategy updates the current supporting set by merging the newly selected vector set. During each iteration step, a bottom-up approach is presented to estimate and update the sparsity factor adaptively and heuristically instead of determining it manually or empirically. In addition, our method incorporates an effective technique for re-evaluating the reliability of all candidates at each iteration of the process, which guarantees the accuracy and reliability for fluorescence reconstruction. To better evaluate the proposed method, we compared it to the IS_L1 method and the StOMP method both in numerical experiments and in vivo mouse experiments [15, 17]. The proposed method is proved to be more accurate, efficient, robust, and reliable for fluorescence reconstruction compared to the contrasting methods, which demonstrates its potential for practical FMT applications.
The outline of the paper is listed as following. In section 2, we present the reconstruction methodology for FMT. In the beginning, the diffusion approximation for the radiative transfer equation is brieñy introduced. Then, the adaptive sparsity subspace pursuit based reconstruction method is elaborately formulated. In section 3, numerical phantom experiments of the proposed method are performed. In section 4, experimental mouse reconstruction in heterogeneous tissue further demonstrates the robustness and feasibility of this proposed method. In section 5, we discuss relative issues and conclude our work.
2.1 Photon propagation model
Describing photon propagation in biological tissues can be modeled using the integro-differential equation known as the radiative transfer equation (RTE) . However, it is a major endeavor to provide solutions for the RTE and it remains a challenging task in the fields of tissue optics and radiological sciences . For photon propagation in biological tissues within the near-infrared spectral window, scattering is the dominant phenomenon over absorption. Therefore, the diffusion equation (DE) as a low-order approximation of the RTE is usually used to model the photon transport in highly scattering media [22, 23]. For steady-state FMT with point excitation sources, the following coupled diffusion equations are widely used as an effective model for photon propagation [23, 24]:
where denotes the domain of the problem; subscripts x and m denote the excitation and emission wavelengths respectively, is the photon flux density for excitation (subscript x) and fluorescent light (subscript m); denotes the absorption coefficient of the tissue and denotes the diffusion coefficient; denotes the scattering coefficient; and g denotes the anisotropy parameter . is the fluorescent yield to be reconstructed which incorporates the quantum efficiency and absorption coefficient of the fluorescent probe. Here, we assume that the absorption and scattering of the excitation light caused by the fluorescent probes will have little effect on the distribution of , because the fluorescent probes usually occupy a very small space compared to domain . represents the different excitation point source positions with the amplitude . For FMT with point illumination, a collimated laser spot is usually modeled as an isotropic point source , where is the point one transport mean free path into the medium from the illumination spot; is the amplitude of the point source; and is the Dirac function.
In order to solve the coupled equations, the Robin-type boundary conditions (RBC) are implemented on the boundary of domain :
where denotes the outward normal vector to the surface . is a constant that is approximated as . is a parameter governing the internal reflection at the boundary .
2.2 Linear relationship establishment
A finite element method is used to calculate solutions of the diffusion Eqs. (1) together with the boundary condition (2). We discretize the domain with small tetrahedrons and take the base functions as the test functions. Then, the FMT problem can be linearized and the following matrix-form equations can be obtained:
where and denote the base function of node i and node j, respectively; denotes the system matrix for excitation (subscript ) and emission (subscript ); and denotes the excitation source distribution after discretization. Matrix is obtained by discretizing the unknown fluorescent yield distribution, while vector X is the fluorescent yield yet to be reconstructed.
For each excitation point source at during the excitation process, the photon density , which is used as the energy source for the emission process, can be obtained by solving Eq. (3). Considering that is symmetrical and positive definite, Eq. (4) can be formulated into the following matrix-form equation:Eq. (6) for different excitation positions and obtain the following Matrix-form equation:4, 7, 26, 27].
Hence the linear relationship between the measured photon density distribution and the unknown fluorescent source distribution has been established, and can improve computational efficiency. It is worth mentioning that Eq. (7) is an underdetermined system of linear equations because the available surface photon density distribution is far more limited than the unknown internal fluorescent source distribution, which causes difficulties in the following tomographic reconstruction.
2.3 Reconstruction based on adaptive sparsity subspace pursuit
The reconstruction for FMT is involved in solving an underdetermined system of linear equations. As mentioned above, because tumors are small and sparse compared to the entire reconstruction domain in earlier detection, the fluorescent sources, which indicate the distribution of tumors, are usually sparse. This can be considered as a kind of a priori information of the fluorescent sources. Here, the L1 regularization is incorporated into the FMT problem to promote the sparsity of the solution, and Eq. (7) is transformed into the one-norm regularized least-squares problem:Eq. (9) is expected to be a sparse solution (i.e., solution where only a small number of entries are nonzero). The SASP method aims to achieve a sparse solution to the FMT problem. Next, we will describe its basic ingredients.
In the SASP method, some innovative strategies including subspace projection, the bottom-up sparsity adaptive approach, and backtracking technique are applied to enhance the reconstruction accuracy and robustness. The main steps of the proposed algorithm are summarized below.
The SASP method is a heuristic method. It starts with an empty index set and initial residual vector . The sparsity factor K and the step size S are also initialized. To further illustrate how such a method fits into FMT reconstruction, we mathematically derive the SASP method using the following five steps.
To make the candidate set of the proposed method by using the maximal correlation test, assuming that within an individual iteration of optimization, we have the (n-1)th estimate for the residual vector . Then, the n-th iteration applies matched filtering to the current residual vector, getting a vector of residual correlations to be:
To form the true support set of the proposed method, we introduce the subspace projection approach and the backtracking technique. Next, the observation vector is projected onto the subspace spanned by the columns of , where denotes the submatrix of with indices from the candidate set . The projection of onto span() is defined as
To determine whether the new support set is well suited to build the sparse solution of the FMT problem, we need to calculate the new residual. We then perform projection of the observations onto the subspace spanned by the columns of belonging to the new support set . Let denote the submatrix with columns chosen using index set , we have
To estimate the optimal sparsity factor of the fluorescent sources, the proposed method adopts a sparsity adaptive strategy, which alternatively estimates the sparsity and the true support set of the fluorescent sources.
Here, an adaptive one-dimensional search approach is used to get the optimal sparsity factor K according to the changes from the previous residual to the new residual . If the norm of is smaller than that of , the current sparsity factor K is considered to be suitable for the reconstruction process in the current iteration and doesn’t need to be changed. The residual and the true support set , which can be used in the next iteration, are replaced by the newly estimated residual and support set respectively. Otherwise, the sparsity factor is estimated to be smaller than the real case, and it is updated by adding the step size , i.e., the new sparsity factor is set as K = K + S. In addition, the residual and true support set are replaced by the previous residual and , ensuring that the newly generated residual is the minimal residual until the current iteration.
This step is used to determine whether or not it is time to discontinue this proposed method. There are many halting conditions for practical FMT reconstruction methods. One common approach is to stop the method when the norm of the current residual is below a certain threshold value. Another approach is to end the method when the relative residue improvement between two consecutive iterations was below a certain threshold, because it is not worth it to take more costly iterations if the resulting improvement is too small. In this paper, we suggest that the SASP method is finished when the norm of the residual is smaller than the threshold or the maximum iteration number is reached. In our following experiments, the parameters and were optimized according to experimental experience. The threshold was set to equal , where denoted the Euclidean length of the observation vector . In our experiments, the proposed method was good enough to obtain satisfactory results when using as the value of the threshold parameter . The maximum iteration number was set to be 25. In our experiments, all the reconstructions with the proposed method stopped within 25 iterations.
We check the halting conditions and, if it is not yet ready to stop, we set and go to the next iteration of the SASP method. In the next iteration, the updated sparsity factor and residual will be employed. If it is time to stop, the current support set is good enough to build the sparse solution of the FMT problem. We then output the estimated sparse solution , satisfying and , where denotes the length of .
To better understand the proposed reconstruction method in our study, further instructions on the backtracking technique are given.
The backtracking technique: Note that there are 2K indices in the candidate set , among which some of them do not belong to the correct support set. In Step II of the SASP method, a backtracking strategy is adopted to expurgate those indices from and form the new support set . This backtracking strategy is quite important as it enables the method to re-evaluate the reliability of all candidates at each iteration of the process, and it can remove the indices that do not belong to the current support set and are added in the previous iteration. Those indices are thought to be right in the previous iteration, but they are thought to be wrong in the current iteration. If some of the removed indices belong to the final support set of the reconstruction method, they will be added to the support set in the following iterations. The backtracking strategy will make the method more robust and make the reconstruction results more accurate. There is a big difference between the StOMP method and the SASP method because the former one generates a list of candidates sequentially, without backtracking. This difference makes the SASP method more reliable than the StOMP method.
3. Numerical experiments and results
In this section, numerical results were presented to demonstrate the feasibility of the proposed reconstruction method. The numerical experiments were performed based on a mouse-mimicking heterogeneous phantom, as shown in Fig. 1. The phantom was a cylinder, which was 20mm in diameter and 20mm in height, consisting of four kinds of materials to represent muscle (M), lungs (L), heart (H), and bone (B). The optical parameters of each kind of materials for both the excitation and emission wavelength are listed in Table 1 [27–29]. Figure 1(a) demonstrates a 3D view of the phantom which was discretized into 3756 nodes and 21715 tetrahedral elements using the finite-element method and Fig. 1(b) demonstrates the cross-section of the phantom in the z = 0 plane. The black dots in Fig. 1(b) represent the excitation light sources, which were modeled as isotropic point sources placed at a depth of one transport mean free path beneath the surface in the z = 0 plane. Fluorescence measurement was conducted in transillumination mode. For each excitation source, the emitted fluorescence was measured from the opposite side of the cylindrical phantom with a 160° field of view (FOV), which is illustrated in Fig. 1(b). As mentioned above, the fluorescent sources are usually quite small and sparse for FMT, so, we used small spheres with a diameter of 2mm centered in the z = 0 plane to represent the fluorescent sources. Figure 2 demonstrates three different phantom setups under the situation of single, double and triple sources respectively. The first row of Fig. 2(a), 2(b), 2(c) shows the 3D views of the phantom setups and the second row of Fig. 2(d), 2(e), 2(f) shows the corresponding slice images in the z = 0 plane. The fluorescent yields of the sources were all set to be 0.6 mm−1. The finite element method was used to solve the forward model and the measurements on the surface of the phantom were collected [27, 29]. Considering that in practical FMT experiments using a CCD camera, the shot noise always exists. When large numbers of photons are collected, the noise will approach a Gaussian distribution. We added 5% Gaussian noise to simulate the real case.
In this paper, in order to better evaluate the proposed method, we compared it with the IS_L1 method and the StOMP method. The second reconstruction method of the two contrast methods was adopted as it used a greedy pursuit strategy. For the proposed method, the step size S was set to 2. For the IS_L1 method, the regularization parameter was set to 2e−5. For the StOMP method, the parameter was set to 0.8 and the parameter was set to 100, which were the same as the values used by D. Han et al. . For all of the reconstruction methods, the zero vector was used as the initial value of the solution. All of the reconstructions were performed on our desktop computer with 2.39 GHz Intel Core 2 Duo CPU and 2GB RAM.
3.1 Experiment 1——evaluation of reconstruction accuracy
In the first experiment, we reconstructed fluorescent sources with different numbers to evaluate the reconstruction accuracy of the proposed method. Three different phantom setups for single, double and triple sources were used, as presented in Fig. 2. The fluorescence was excited by point sources from 12 different locations in sequence, as presented in Fig. 1(b). Every 30°, measurements of emission fluorescence were taken in the transillumination mode, then, a total number of 12 data sets were collected to perform the reconstruction of the fluorescent yield. In order to simulate the real case, 5% Gaussian noise was added to the measurement data sets. The actual fluorescence yield was set as 0.6 mm−1. To validate the reconstruction accuracy of the proposed method, we compared it to the IS_L1 method and StOMP method. Figures 3, 4 and 5 give the reconstruction results for three different phantom setups, demonstrating the 3D view of the reconstructed sources in the phantom setups and the slice images in the z = 0 plane. The red circles in the slice images denote the real locations of the fluorescent sources. As shown in Figs. 3, 4 and 5, the stable locations of the fluorescent sources can be obtained through IS_L1 method, but the maximum values of the reconstructed fluorescent yield were smaller than those of the StOMP method and the SASP method. The StOMP method could acquire acceptable results when only one or two fluorescent sources existed. However, when three fluorescent sources existed, the reconstruction results were not satisfactory. Based on the data in Fig. 5(b) and Fig. 5(e), the fluorescent yield of the two sources on the left side (i.e., Source S1 and S2) was not well reconstructed. The reason for this result may be that the StOMP method didn’t adopt a sparsity adaptive strategy to achieve the optimal sparsity factor of the fluorescent sources and the choices of threshold parameters and were estimated before the reconstruction process started and mainly came from experience. The proposed method could obtain satisfactory reconstruction results in three cases. These more accurate results benefited from the backtracking strategy and sparsity adaptive strategy. To analyze the reconstruction accuracy quantitatively, we defined the relative intensity error to be , where denotes the true fluorescent yield of the source and denotes the maximum reconstructed value of the fluorescent yield. We also defined the location error to be , where denotes the true location of the source center and denotes the location of the finite element node with the maximum reconstructed value of the fluorescent yield for that source. Table 2 summarizes the results of the quantitative analysis of different reconstruction methods. As shown in Table 2, all the three reconstruction methods could acquire satisfactory source localizations; however, the relative intensity errors of the proposed method were smaller. Figure 6(a) shows the relationship curve between sparsity factor K and iteration step n with a different number of sources. The sparsity factor K was updated automatically according to the changes of the residual between two consecutive iterations. Simultaneously, the residual was reduced to a small value, as shown in Fig. 6(b). If the residual was below , the relative residual decrement between two consecutive iterations became too small and we thought it was not worth to take more costly iterations.
3.2 Experiment 2——evaluation of reconstruction efficiency
To examine the reconstruction efficiency of the proposed method, we also compared it with the IS_L1 method and the StOMP method, which were two mainstream methods for FMT reconstruction. This numerical experiment was performed using a mouse-mimicking heterogeneous cylindrical phantom which contained three sources, as presented in Fig. 2 (c, f). To better evaluate the time cost of these reconstruction methods, we conducted the numerical experiment on five kinds of volumetric meshes of different sizes. For all these reconstruction methods, the zero vector was used to initialize the unknowns. The comparison of the reconstruction efficiency was made based on the data in Table 3. It showed the time consumed by the proposed method as well as the two other methods compared for reconstructing the five groups of data sets whose sizes were determined by the density of the discrete volumetric mesh. Each value of the reconstruction time was the average of ten independent runs. All three reconstruction methods were coded in MATLAB and performed on our desktop computer with 2.39 GHz Intel Core 2 Duo CPU and 2GB RAM.
The results revealed that: 1. When reconstructing the same data set, the StOMP method worked more efficiently than the IS_L1 method, but was less efficient than the proposed SASP method; 2. When the size of the data set was increased, the SASP method became more computationally competitive; 3. All of the data sets were discretized based on the mouse-mimicking heterogeneous phantom with three sources in the tissue material of the lungs (L) (Fig. 2 (c, f)), which implied that the proposed method was not only time efficient but also potential for practical FMT applications.
3.3 Experiment 3——evaluation of reconstruction robustness of limited measurement data
To verify the reconstruction robustness of the proposed method, the amount of measurement data was reduced and the reconstruction of the fluorescent sources was conducted, simulating a much worse case scenario. This is quite necessary because long-term measurement is not appropriate or feasible in some cases of practical FMT experiments. For example, when imaging small animals such as mice or rabbits, the bleaching effect of the fluorescent probe caused by long-term measurement must be taken into consideration, as it can affect the experimental results. To resolve this problem, it is necessary to reduce the amount of fluorescence measurements and reconstruct the fluorescent sources from quite limited measurement data. Here, only the measurement data sets generated by excitation point sources 1, 4, 7 and 10 were retained by us, as shown in Fig. 1(b). In other words, we conducted the reconstruction from only four measurement data sets corrupted by 5% Gaussian noise. A mouse-mimicking heterogeneous cylindrical phantom with three sources was used to perform the numerical experiment, as shown in Fig. 2(c), 2(f). Figure 7 presents the reconstruction results from four measurement data sets, demonstrating the 3D views of the reconstructed sources in the phantoms and the corresponding slice images in the z = 0 plane. The location errors of the reconstruction results, as well as the relative intensity errors, are summarized in Table 4.
As shown in Fig. 7 and Table 4, the reconstructed sources by the StOMP method were not accurately located in the true regions. The IS_L1 method could obtain satisfactory source localization for the fluorescent sources S1 and S2, but the location error for the reconstructed source S2 was large. On the contrary, the proposed method could recover all three fluorescent sources accurately. In addition, the relative intensity errors for the proposed method were smaller than the other two compared methods. The above results suggested that the proposed method could obtain satisfactory results even under quite ill-posed conditions, and had the potential for practical FMT applications.
4. In vivo experiments
In this section, to further study the potential of the proposed method in the practical application of FMT, an in vivo experiment on an adult Kunming mouse (Laboratory Animal Center, Peking University, China), which was implanted with a plastic fluorescent bead in the muscle has been conducted using the hybrid optical/micro-CT imaging system developed by our group [30–32]. This hybrid imaging system, surrounded in a completely dark environment, provides two modalities: multi-view FMT and micro-CT. The schematic illustration of the system is shown in Fig. 8, which is mainly equipped with a cooled Charge Coupled Device (CCD) camera (Princeton Instruments, PIXIS: 1024), a 671mn continuous wave (CW) laser, a set of optical lenses, an x-ray generator (Oxford Instruments, 90 kV UltraBright Micro-focus Source), an x-ray detector (Hamamatsu Photonics, C7943CA-02 Flat Panel Sensor), and a rotating stage (Beijing Zolix Instruments, RAK). The CCD camera is used as an optical detector which was cooled to −70°C. The excitation light source was provided by a 671 nm CW laser with an output power of 22 mW and a laser spot diameter of 1mm, which simulated a point source case. The fluorescence measurements were collected in transillumination mode. A band-pass filter with a bandwidth of 10 nm and a center wavelength of 700nm was used to allow light transmission at the emission wave-length.
The main process of in vivo experiments can be summarized as following. Firstly, the mouse was injected through the caudal tail vein with 0.4 ml Fenestra vascular contrast agent to enhance the micro-CT images, followed by 0.3 ml of anesthetic at a 0.15g/ml concentration via intraperitoneal injection. When the anesthetic set in, a bead filled with cy5.5 solution with a concentration of 2000nM was embedded stereotactically into the body of the mouse in the vicinity of its liver. This fluorescent solution had an extinction coefficient of 0.019 and a quantum efficiency of 0.23 at the peak excitation wavelength of 671nm . It is worth to mention that the reconstruction accuracy of the proposed method can be examined in this study, since the luminescent bead was wrapped in a plastic material, which could be easily detected by the micro-CT system. About half an hour after the injection of the Fenestra vascular contrast agent, the mouse was placed on the automatic rotation stage, as delineated in Fig. 8. Secondly, we began to collect the dual-modality raw data respectively. The acquisition time for the entire procedure was within 10min. The fluorescent signals emitted from the surface of the mouse body were captured in several different views by the cooled CCD camera. And then the cone-beam micro-CT projection data with 360° views was scanned using the x-ray generator and detector.
After collecting the raw data, some essential preprocessing operations were carried out and set the data set suitable for posterior FMT reconstruction.
- (1). Micro-CT Reconstruction: The cone-beam micro-CT projection data were reconstructed by the Feldkamp-Davis-Kress (FDK) algorithm on commodity GPU using an acceleration scheme  to yield the 3D volume data, three slices of which are shown in Fig. 9 where the yellow squares mark the location of the fluorescent bead at the coordinate (34.40, 33.80, 6.40)
- (2). Segmentation: In order to build the heterogeneous mouse model delineated in Fig. 10, we performed segmentation on six main kinds of organs including muscle, lungs, heart, liver, bone, and kidneys through a combination method of interactive and threshold segmentation. The optical properties for different organs were calculated according to the literature  and are listed in Table 5
- (3). Discretization: Prior to solving the FMT inverse problem using the finite element method, we performed meshing on the segmented micro-CT data of the mouse body. Here, the torso section of the mouse was discretized into a volumetric mesh containing 4,659 nodes and 22,667 tetrahedral elements
- (4). Registration: In order to portray the photon distribution on the surface of the mouse body, the optical fluorescent images were mapped onto the surface of the volumetric mesh in terms of space and energy. The surface energy mapping was carried out by using a 3D surface flux reconstruction algorithm , as is seen in Fig. 10(c)
After the above procedures, we used the finite element method to solve the FMT inverse problem and used three different reconstruction methods (the IS_L1 method, the StOMP method, and the proposed method) to perform the FMT reconstruction. The 3D results reconstructed by the IS_L1 method, the StOMP method and the proposed method are shown in Fig. 11. The muscle region is set to be translucent, so that the reconstructed sources is not covered. Quantitative comparisons of the reconstruction results for the above three methods are shown in Fig. 12 and Table 6. It can be perceived that there were more than one reconstructed sources inside the torso of mouse reconstructed by the StOMP method, and the smaller reconstructed sources were far away from the center of the real source, although the reconstructed location center was (34.05, 33.50, 6.07) and had a location error of only 0.57mm. The reconstructed center position of the fluorescent source was (34.11, 126.96.36.199) for the IS_L1 method, and (34.05, 33.50, 6.07) for the proposed method. The corresponding location errors were 0.57 mm and 1.56 mm respectively. Although both of the reconstructed fluorescent sources had a small distance from the actual source in center position, artifacts appeared near the reconstructed fluorescent source in the results of the IS_L1 method, whereas artifacts did not appear in the proposed method’s results. The reconstructed location error by the proposed method was acceptable, because the errors may have been caused by the diffusion equation model, the geometrical approximation, the optical properties’ inaccuracy, the surface energy mapping, and so on. The reconstruction time for the proposed method was 5.5625s and was much faster than the two contrasting methods, which demonstrated the absolute advantages of the proposed method in efficiency.
The above findings suggested that the proposed method was able to reconstruct the fluorescent source accurately and had the potential to detect the lesions in practical biomedical applications.
5. Discussion and conclusion
In this paper, a novel sparsity adaptive subspace pursuit method is proposed and analyzed for FMT reconstruction. This reconstruction method is proven to be accurate, robust and high-efficient for fluorescent source reconstruction mainly through using the following three strategies:
- (1) A subspace projection and correlation maximization approach has been utilized to simplify the FMT problem and treat the FMT problem with sparsity-promoting L1-norm regularization as the basis pursuit problem. Based on this, a new sparse reconstruction problem, with a greatly reduced order compared to the original optimization problem is formed between the measurement data sets and the permissible coefficients.
- (2) Instead of using traditional methods, which require the sparsity factor as a priori for exact reconstruction, a sparsity adaptive strategy has been proposed to estimate the sparsity of the FMT problem. This innovative feature makes the proposed method a promising candidate for many practical FMT applications when the number of non-zero coefficients of the fluorescent source distribution is not available.
- (3) A backtracking technique has been proposed to re-evaluate the reliability of all candidates at each iteration of the reconstruction process and to remove the wrong items added in the previous iteration, which can make the reconstruction method more robust and the reconstruction results more accurate.
Moreover, three numerical experiments based on a mouse-mimicking heterogeneous phantom and an in vivo experiment based on an adult Kunming mouse have been conducted to evaluate the feasibility of our method in fluorescent source reconstruction. The results demonstrate that:
- (1) The reconstruction accuracy of the proposed method has been validated by the experiment on different mouse-mimicking heterogeneous phantom setups for a single source, double sources, and triple sources, and the results indicate that it is able to localize different fluorescent sources with a position bias less than 1 mm.
- (2) The proposed reconstruction method has been proven to be more efficient compared to the classical approaches, and it becomes more computationally competitive as the data set grows larger.
- (3) The reconstruction robustness of the proposed method has been verified in that it can obtain satisfactory fluorescent source reconstruction results even under quite ill-posed conditions, for example, when the measurement data sets are quite limited.
- (4) The potential of the proposed method on the application of FMT has been further validated by an in vivo experiment on an adult Kunming mouse model, in which a small fluorescent source in the vicinity of the mouse’s liver has been accurately reconstructed.
For FMT problems, the diffuse approximation model has been extensively used to mimick photon propagation in biological tissues. The advantage of this model is computationally efficient and has an explicit physical meaning. However, more accurate models to describe photon propagation in biological tissues are also necessary for higher quality of 3D imaging reconstruction. In recent years, although the computational complexity of these models are high, several complex models such as the simplified spherical harmonics approximation (SPN) model have been proposed to improve the optical tomography reconstruction quality. Our method can also be utilized in these complex models with few modifications.
Future work will focus on in vivo experiments with probe-marked mouse models, which is quite important for tumor detection. In addition, the in vivo experiment will be extended further to detect weak optical signals from the internal fluorescent sources at different depths as well as multiple sources instead of only one. The FMT reconstruction for mouse models with multiple sources will be conducted in future studies.
This paper is supported by the National Basic Research Program of China (973 Program) under Grant No. 2011CB707700, the National Natural Science Foundation of China under Grant Nos. 81027002 and 81227901, and the Chinese Academy of Sciences Visiting Professorship for Senior International Scientists under Grant No. 2012T1G0036. The authors thank Adjunct Assistant Professor Dr. Karen M. von Deneen from the University of Florida and Dr. Yang Du for their modifications in this paper.
References and links
1. H. Fan-Minogue, Z. Cao, R. Paulmurugan, C. T. Chan, T. F. Massoud, D. W. Felsher, and S. S. Gambhir, “Noninvasive molecular imaging of c-Myc activation in living mice,” Proc. Natl. Acad. Sci. U.S.A. 107(36), 15892–15897 (2010). [CrossRef] [PubMed]
2. M. A. Whitney, J. L. Crisp, L. T. Nguyen, B. Friedman, L. A. Gross, P. Steinbach, R. Y. Tsien, and Q. T. Nguyen, “Fluorescent peptides highlight peripheral nerves during surgery in mice,” Nat. Biotechnol. 29(4), 352–356 (2011). [CrossRef] [PubMed]
3. G. M. van Dam, G. Themelis, L. M. Crane, N. J. Harlaar, R. G. Pleijhuis, W. Kelder, A. Sarantopoulos, J. S. de Jong, H. J. Arts, A. G. van der Zee, J. Bart, P. S. Low, and V. Ntziachristos, “Intraoperative tumor-specific fluorescence imaging in ovarian cancer by folate receptor-α targeting: first in-human results,” Nat. Med. 17(10), 1315–1319 (2011). [CrossRef] [PubMed]
5. A. Ale, V. Ermolayev, E. Herzog, C. Cohrs, M. H. de Angelis, and V. Ntziachristos, “FMT-XCT: in vivo animal studies with hybrid fluorescence molecular tomography-X-ray computed tomography,” Nat. Methods 9(6), 615–620 (2012). [CrossRef] [PubMed]
7. X. L. Song, D. F. Wang, N. G. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007). [CrossRef] [PubMed]
9. W. Bangerth and A. Joshi, “Adaptive finite element methods for the solution of inverse problems in optical tomography,” Inverse Probl. 24(3), 034011 (2008). [CrossRef]
10. N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007). [CrossRef] [PubMed]
11. Y. Lin, H. Gao, O. Nalcioglu, and G. Gulsen, “Fluorescence diffuse optical tomography with functional and anatomical a priori information: feasibility study,” Phys. Med. Biol. 52(18), 5569–5585 (2007). [CrossRef] [PubMed]
12. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]
13. P. Mohajerani, A. A. Eftekhar, J. D. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46(10), 1679–1685 (2007). [CrossRef] [PubMed]
14. D. Han, J. Tian, K. Liu, J. Feng, B. Zhang, X. Ma, and C. Qin, “Sparsity-promoting tomographic fluorescence imaging with simplified spherical harmonics approximation,” IEEE Trans. Biomed. Eng. 57(10), 2564–2567 (2010). [CrossRef] [PubMed]
15. D. Han, J. Tian, S. P. Zhu, J. C. Feng, C. H. Qin, B. Zhang, and X. Yang, “A fast reconstruction algorithm for fluorescence molecular tomography with sparsity regularization,” Opt. Express 18(8), 8630–8646 (2010). [CrossRef] [PubMed]
16. M. Elad, B. Matalon, and M. Zibulevsky, “Coordinate and subspace optimization methods for linear least squares with non-quadratic regularization,” Appl. Comput. Harmon. Anal. 23(3), 346–367 (2007). [CrossRef]
17. D. Han, X. Yang, K. Liu, C. Qin, B. Zhang, X. Ma, and J. Tian, “Efficient reconstruction method for L1 regularization in fluorescence molecular tomography,” Appl. Opt. 49(36), 6930–6937 (2010). [CrossRef] [PubMed]
18. N. Blow, “In vivo molecular imaging: the inside job,” Nat. Methods 6(6), 465–469 (2009). [CrossRef]
19. A. M. Bruckstein, D. L. Donoho, and M. Elad, “From sparse solutions of systems of equations to sparse modeling of signals and images,” SIAM Rev. 51(1), 34–81 (2009). [CrossRef]
20. A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. 202(1), 323–345 (2005). [CrossRef]
21. A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. 220(1), 441–470 (2006). [CrossRef]
24. J. H. Lee, A. Joshi, and E. M. Sevick-Muraca, “Fully adaptive finite element based tomography using tetrahedral dual-meshing for fluorescence enhanced optical imaging in tissue,” Opt. Express 15(11), 6955–6975 (2007). [CrossRef] [PubMed]
25. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The Finite Element Method for the Propagation of Light in Scattering Media: Boundary and Source Conditions,” Med. Phys. 22(11), 1779–1792 (1995). [CrossRef] [PubMed]
26. D. F. Wang, X. Liu, Y. P. Chen, and J. Bai, “A Novel Finite-Element-Based Algorithm for Fluorescence Molecular Tomography of Heterogeneous Media,” IEEE Trans. Inf. Technol. Biomed. 13(5), 766–773 (2009). [CrossRef] [PubMed]
28. A. B. Milstein, J. J. Stott, S. Oh, D. A. Boas, R. P. Millane, C. A. Bouman, and K. J. Webb, “Fluorescence optical diffusion tomography using multiple-frequency data,” J. Opt. Soc. Am. A 21(6), 1035–1049 (2004). [CrossRef] [PubMed]
29. W. X. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005). [CrossRef] [PubMed]
32. P. Wu, K. Liu, Q. Zhang, Z. Xue, Y. Li, N. Ning, X. Yang, X. Li, and J. Tian, “Detection of mouse liver cancer via a parallel iterative shrinkage method in hybrid optical/microcomputed tomography imaging,” J. Biomed. Opt. 17(12), 126012 (2012). [CrossRef] [PubMed]
33. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography,” Opt. Express 14(16), 7109–7124 (2006). [CrossRef] [PubMed]
34. G. Yan, J. Tian, S. Zhu, Y. Dai, and C. Qin, “Fast cone-beam CT image reconstruction using GPU hardware,” J. XRay Sci. Technol. 16, 225–234 (2008).
35. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [CrossRef] [PubMed]
36. X. Chen, X. Gao, D. Chen, X. Ma, X. Zhao, M. Shen, X. Li, X. Qu, J. Liang, J. Ripoll, and J. Tian, “3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images,” Opt. Express 18(19), 19876–19893 (2010). [CrossRef] [PubMed]