Hybrid Multilevel Sparse Reconstruction for a Whole Domain Bioluminescence Tomography Using Adaptive Finite Element

Quantitative reconstruction of bioluminescent sources from boundary measurements is a challenging ill-posed inverse problem owing to the high degree of absorption and scattering of light through tissue. We present a hybrid multilevel reconstruction scheme by combining the ability of sparse regularization with the advantage of adaptive finite element method. In view of the characteristics of different discretization levels, two different inversion algorithms are employed on the initial coarse mesh and the succeeding ones to strike a balance between stability and efficiency. Numerical experiment results with a digital mouse model demonstrate that the proposed scheme can accurately localize and quantify source distribution while maintaining reconstruction stability and computational economy. The effectiveness of this hybrid reconstruction scheme is further confirmed with in vivo experiments.


Introduction
Bioluminescence imaging (BLI) is an in vivo imaging modality that has been successfully used in preclinical researches [1][2][3]. This imaging strategy exploits the properties of luciferase that can generate visible or near infrared light through the oxidation of an enzyme-specific substrate in the presence of oxygen and adenosine triphosphate [4]. As the produced light intensity is directly proportional to the concentration of luciferase-expressing cells, BLI can reveal cellular and molecular features of biology and disease [5]. However, BLI fails to provide depth information of the internal biological sources [6]. Collecting measurement data from multiple views or combining multiple BLI acquisition with geometrical structures acquired by micro-CT or MRI, bioluminescence tomography (BLT) tries to reconstruct the 3D biological source distribution. In this way, BLT overcomes the limitation of planar imaging in poor spatial resolution and further facilitates our understanding of biomolecular processes as they occur in living animals. Therefore, BLT has substantial potential to be a powerful tool for noninvasively monitoring and tracking a variety of biological processes [7].
Generally, BLT involves a forward and an inverse problem (source reconstruction). Due to the diffusive nature of photon propagation in tissue, BLT source reconstruction is known to be a highly ill-posed problem [6,8]. To overcome the inherent ill-posedness of the tomographic problem in BLT, different strategies have been proposed either by increasing the amount of independent measurements with spectrally resolved or multispectral approaches [9][10][11][12][13] or by reducing the number of unknowns with permissible source region [6,10,13,14]. Up to now, quantitative reconstruction for whole domain BLT with monochromatic boundary measurements has not been intensively investigated.
As in many other imaging modalities, the achievable resolution for BLT is determined firstly by the signal to noise ratio, and secondly by the level of discretization. Image quality can be improved by uniformly refining mesh throughout the reconstruction domain. Nevertheless, global refinement tends to further aggravate the ill-posedness and incur insurmountable computational burden due to the increased unknowns and problem size. Consequently, the use of adaptive finite element method (AFEM) is an indispensable approach to improve image quality [15][16][17][18][19][20][21].
In this contribution, we present a whole domain BLT method based on AFEM which provides fine resolution around targets with coarser resolution in other regions. Unlike the previous AFEM-based BLT that adopted identical inversion strategy on different mesh levels [15,[18][19][20][21], we take the variance on different discretization levels into account and propose a novel hybrid multilevel reconstruction scheme to maintain solution stability and computational economy. Two different inversion algorithms, the stagewise fast LASSO (SwF-LASSO) [22] and the incomplete variables truncated conjugate gradient method (IVTCG) [23], are applied to the first mesh level and the succeeding ones according to their respective characteristics.
The following sections describe some of the implementation details of the hybrid AFEM algorithm, the evaluations on a digital mouse model, and the validation with an in vivo experiment. Short discussions and concluding remarks are given at the end of this paper.

Photon Propagation Model.
In this work, we assume that the structural and optical parameters regarding different organs are given. Therefore, the BLT reconstruction comes down to a linear inverse source problem. Based on the diffusion approximation model of radiative transfer equation, a linear relationship between the source distribution and boundary measurements is then derived with the finite element method [6]: where ∈ × ( < ) is the system matrix, ∈ denotes the internal source distribution, and Φ ∈ represents measurable boundary nodal photon density that is usually calculated from the surface flux image captured by a CCD camera.
In view of the limitation of using permissible source region in BLT reconstruction, we consider a whole domain reconstruction scheme without this kind of a priori information. On the other hand, 1 -norm based sparse regularization methods have attracted considerable amount of attention in BLT [10,[20][21][22][23][24][25], and the reconstructions' results therein demonstrate that 1 -norm solution fits the sparsity nature of bioluminescent source distribution in BLT practice. Using 1 regularization, we formulate the BLT inverse problem to the following optimization problem: where ‖ ⋅ ‖ 2 denotes the Euclidean norm, ‖ ⋅ ‖ 1 is the 1 norm, and > 0 is a regularization parameter.

Hybrid Multilevel Reconstruction Based on AFEM.
In order to provide the resolution necessary for imaging at acceptable computational cost, the domain Ω is dynamically discretized into a nested sequence of tetrahedral meshes {Θ 1 , . . . Θ , . . .}, rather than a fixed and uniformly fine mesh.
In the proposed hybrid multilevel AFEM reconstruction process, reconstruction starts at the coarsest level and proceeds to the finer ones by locally refining the particular region based on a previous reconstruction procedure. We note that the first reconstructed procedure on the coarsest mesh is quite different from the subsequent ones in the following three aspects. (i) It is based on a uniform mesh while others are with a locally refined mesh. (ii) The inversion on the first discretization level involves a large-size underdetermined system. In contrast, all of the subsequent reconstructions on locally finer region involve overdetermined systems. (iii) It has no a priori information of a promising region in whole domain case, whereas the others can obtain a permissible source region to constrain the solution space from a previous reconstruction procedure. Consequently, the specific inversion should be different on different meshes, and thus we propose a hybrid multilevel reconstruction scheme.
On the first mesh Θ 1 , we employed the recently reported greedy algorithm SwF-LASSO to solve the underdetermined problem in (2). The SwF-LASSO algorithm converges very fast and is able to find an approximate value close to the real distribution in only a few iteration steps. A brief outline of SwF-LASSO is given as follows [22].
where = , = ( − ) −1 and consists of those column vectors of A relating to the selected basis functions in +1 . The updating formula of S is where Step 4. = − +1 , and = + +1 .
After the inversion on Θ ( = 1, 2, 3, . . .) completes, adaptive mesh refinement is triggered. All of the elements Obtain the next locally finer mesh Θ with the adaptive mesh refinement strategy.
Take all of the nodes with nonzero value into account to form a permissible source region and build a reduced model: Solve the model with the IVTCG method and obtain the reconstruction result on this locally refined mesh Θ  with nonzero reconstructed value are selected to be refined, which can be regarded as a kind of mesh refinement strategy based on posteriori error estimation. Using the longest-edge bisection method, a locally refined mesh Θ +1 is obtained [17]. Unlike other previous reports, we employ a different reconstruction procedure on the succeeding mesh levels Θ ( > 1). The IVTCG algorithm proposed in [23] has been demonstrated as an effective reconstruction method by reformulating (2) as a convex quadratic program with nonnegative constrained conditions. It updates only partial variables in working set per iteration and adopts a working set splitting strategy to find the searching direction more efficiently, which leads to a small subproblem to be minimized and greatly decreases the number of iterations. The model transformation and the mechanism of IVTCG are detailed in [23].
We note that it is the sparseness-related parameter that controls the size of the subproblem, which is solved by the truncated conjugate gradient method. Generally, for a very sparse problem, IVTCG can obtain accurate results with reasonable computational efficiency by setting = ⌊ /10⌋ and the maximum iterate number of the subproblem iter max = . However, in the reconstruction procedures after local mesh refinement, the target is not a very sparse signal and the computational cost will increase sharply. In view of this feature, we make a modification and adjust the parameter = ⌊ /4⌋, and iter max = 25 in our implementation.
A new round of local mesh refinement and reconstruction will be performed until the number of refinement exceeds the maximum number max or the model misfit ‖ − Φ ‖ 2 2 is reduced below a prespecified threshold . For the results reported in this work, we used max = 4 and = 10 −5 .
The procedure of the proposed hybrid multilevel reconstructions method is illustrated in Figure 1

Numerical Experiments and Results
We tested the proposed hybrid multilevel reconstruction method with a digital mouse model employing synthetically generated data. In the following simulations, we employed a 3D mouse atlas of CT and cryosection data to provide anatomical information [26]. The CT slices of the mouse were segmented into major anatomical components, including lungs, a heart, a liver, a stomach, kidneys, and muscles. The corresponding optical properties were the same as the settings in [27], as shown in Table 1. The whole region included the mouse torso with a height of 45 mm. In the following numerical experiments, the torso model was discretized into a tetrahedral-element mesh, and synthetic measurements were generated by solving the forward model with FEM. To simulate the noise involved in real BLT experiments, 15% Gaussian white noise was added to the synthetic data. The qualities of the reconstruction are quantitatively assessed in terms of location error (LE) and relative error (RE) between the reconstructed power and the actual value.

Single-Target Reconstruction.
In the first set of experiments, a cylindrical source with 0.4 mm radius and 1 mm height was positioned in the right kidney with the center at (11,6,25), as shown in Figure 2(a). The actual source power was 0.2299 nW after discretization with FEM. Figure 2(b) shows the initial mesh for reconstruction and the photon distribution on the surface.
Following the proposed hybrid multilevel reconstructions method, the final result in single-source case was obtained by four rounds of reconstructions. Figure 3 shows the refinement of local mesh around targets and the solution progress from mesh Θ 1 to mesh Θ 4 . According to the proposed methods, fine resolution only presents around targets, while coarser resolution retains in other regions, which contributes to reaching the desirable resolution at acceptable computational cost. Figures 4 and 5 are the transverse views and 3D views of reconstruction results from mesh Θ 1 to mesh Θ 4 , which illustrate the improvement of results during adaptive mesh refinement. To demonstrate the necessity and effectiveness of the hybrid reconstruction scheme, we first compared the SwF-LASSO and IVTCG method on the initial coarse mesh Θ 1 , and then we compared the results of hybrid method, that is, SwF-LASSO + IVTCG, with that of only using SwF-LASSO on the succeeding mesh levels. The detailed reconstruction results are presented in Table 2. Obviously, the reconstruction results by IVTCG are inferior to that of SwF-LASSO on Θ 1 , and hybrid AFEM scheme performs better than the traditional AFEM that uses monoalgorithm of SwF-LASSO on the subsequent mesh Θ 2 to mesh Θ 4 .
Owing to the hybrid multilevel reconstruction scheme, the location error and the relative error of power distinctly decrease with the adaptively local mesh refinement. Especially, significant improvement of reconstructed density and power can be seen from the results in Table 2 and Figure 4.

Double-Source Reconstruction.
We also investigated the resolving ability of the proposed method with two closely separated sources. Two cylindrical sources, same as that in the above single-target setting, were located in the right kidney with their centers at (9, 6.5, 25) and (12,4,25), respectively. They were identical in size and density, but the initial powers of them were 0.2120 nW and 0.2250 nW, mainly due to the influence of the mesh. The source setting and the simulated photon distribution are shown in Figure 6.
In double-source case, the multilevel reconstruction terminated on the third mesh level Θ 3 . Figure 7 Table 3. Figure 7 witnesses an apparent advantage of using hybrid scheme. Although the first-round result was biased towards a node between the targets on mesh Θ 1 , the proposed method successfully identified the two targets finally, which should be attributed to both the AFEM and the hybrid strategy. By contrasting Figure 7(d) with Figure 7(e), we can observe that the improvement caused by multilevel reconstruction with AFEM is evident. Nevertheless the final result of using monoalgorithm of SwF-LASSO is obviously inferior to that of using hybrid algorithm in terms of location accuracy and the reconstructed power error. Take source 1 for instance,  the LE of the hybrid scheme reduces by 0.19 mm and the RE of power falls down to 2.3%. As for source 2, the proposed hybrid reconstruction method yields a 78% plunge in relative error of power.

In Vivo Experiments
To further validate the proposed method, an in vivo experiment was performed on an adult nude mouse. The animal procedures were in accordance with the Fourth Military Medical University that approved the animal protocol.
In the in vivo experiment, a capillary approximately 1.25 mm in diameter and 4.08 mm in length was inserted into the abdomen of the nude mouse. The capillary filled with 5 L luminescent liquid served as the testing source in this experiment. The luminescent solution was extracted from a red luminescent light stick (Glow products, Victoria, Canada), and the generated luminescent light had an emission peak wavelength of about 644 nm. The initial total power was 300 nW (the total power = luminescent solution volume × luminescent solution flux density = 5 L × 60 nW/ L). This set of BLT experiments were conducted with a dual-modality BLT/micro-CT system [23]. The anesthetized mouse was first photographed, and luminescent images were taken by a calibrated CCD camera from four directions at 90 degree intervals with different exposure times. The multiview superimposed photographs and luminescent images are shown in Figures 8(a)-8(d).
After the optical data were acquired, the intact mouse was scanned using the Micro-CT. Because of the limited field of view, only the torso section was scanned. The volume data were reconstructed using GPU-accelerated FDK algorithm [28]. From the CT slices, we located the center coordinate (21.44, 27.52, 9.76) of the actual luminescent source. The mouse body was segmented into five anatomical components, including muscle, heart, lungs, liver, and kidneys. The relevant optical properties of the mouse are listed in Table 4 [29].
Based on the collected multiview luminescent images and the volume data of CT, the 3D surface distribution is determined by the mapping algorithm described in [30], as shown in Figure 8(e). After the mapping process, three rounds of  reconstructions on gradually refined meshes were performed with the proposed hybrid method. The reconstruction result on mesh Θ 1 is presented in Figure 9, where the source center is (20.39, 27.98, 9.78) with a deviation of 1.15 mm to the actual center. From mesh Θ 1 to mesh Θ 3 , the source locations are identical, which means that the SwF-LASSO algorithm yields relatively accurate location from the begging. However, the preliminary reconstruction on the initial coarse mesh Θ 1 possesses relative bigger errors in source power. After two rounds of local mesh refinement, the final results of hybrid method  improved prominently. Specifically, the reconstructed power increased from 149.01 nW to 214.60 nW, and the RE of power decreased from 50.33% to 28.47%. The 3D views of the corresponding results on mesh Θ 1 to mesh Θ 3 are presented in Figure 10.

Discussion and Conclusion
We present a novel multilevel reconstruction method for whole domain BLT, which combines the merit of sparse regularization with the advantage of adaptive FEM. Numerical experiment results employing synthetic data with a digital mouse model illustrate that the proposed hybrid multilevel reconstruction scheme is able to accurately localize and quantify source distribution without a priori information of permissible source region and multispectral measurements. The in vivo experiments conducted on a nude mouse with a dual-modality BLT/micro-CT system further validate the proposed method. From the above experiments, we can find that the inversion algorithm on the initial coarse mesh has more important impact on the final result in the proposed hybrid scheme. The SwF-LASSO algorithm is able to provide a good initial localization with better numerical stability, which guides the subsequent reconstruction on finer meshes to obtain more accurate location and power. Furthermore, the experimental results also demonstrate that the hybrid strategy works. Compared with the multilevel reconstruction using monoalgorithm, the hybrid scheme performs better especially for multiple targets reconstruction. Therefore, it is also possible to form another qualified hybrid scheme using some other promising inversion algorithms.
For the sake of computational efficiency, the reconstructions presented in this paper are based on the diffusion equation. Therefore, the inadequately accurate forward model also leads to some inevitable error. The reconstruction performance might be further improved by using more accurate models, which is also the direction of our further work.
In addition to the many advantages of adaptive finite element methods, such as providing fine resolution around targets with coarser resolution in other region, the proposed hybrid scheme has two remarkable features. (i) Reconstruction result evolves adaptively with iterations, and the reconstruction accuracy is easily controlled by users. (ii) The inversion techniques employed on the initial coarse mesh and the succeeding ones vary with the discretization level to maintain solution stability and computational efficiency.