Multi-atlas registration and adaptive hexahedral voxel discretization for fast bioluminescence tomography

: Bioluminescence tomography (BLT) has been a valuable optical molecular imaging technique to non-invasively depict the cellular and molecular processes in living animals with high sensitivity and specificity. Due to the inherent ill-posedness of BLT, a priori information of anatomical structure is usually incorporated into the reconstruction. The structural information is usually provided by computed tomography (CT) or magnetic resonance imaging (MRI). In order to obtain better quantitative results, BLT reconstruction with heterogeneous tissues needs to segment the internal organs and discretize them into meshes with the finite element method (FEM). It is time-consuming and difficult to handle the segmentation and discretization problems. In this paper, we present a fast reconstruction method for BLT based on multi-atlas registration and adaptive voxel discretization to relieve the complicated data processing procedure involved in the hybrid BLT/CT system. A multi-atlas registration method is first adopted to estimate the internal organ distribution of the imaged animal. Then, the animal volume is adaptively discretized into hexahedral voxels, which are fed into FEM for the following BLT reconstruction. The proposed method is validated in both numerical simulation and an in vivo study. The results demonstrate that the proposed method can reconstruct the bioluminescence source efficiently with satisfactory accuracy.


Introduction
Bioluminescence tomography (BLT) is a non-invasive optical imaging technique used to depict the molecular, cellular and genetic processes in vivo [1][2][3]. Compared with other imaging techniques, BLT has the advantages of low cost, high sensitivity, and high specificity. Therefore, BLT has been widely used in biomedical imaging applications such as monitoring tumor growth or tracking labeled cells [4,5].
BLT is capable of three dimensional (3D) recovery of the location and concentration of the optical probe inside a small living animal from the detected light signal on the surface. However, BLT is a severely undetermined and inherent ill-posed problem, which usually requires additional information to obtain meaningful quantitative reconstruction [6][7][8][9][10]. The optical property uncertainty of the tissue also has a great effect on BLT reconstruction [11]. The structural information provided by computed tomography (CT) or magnetic resonance imaging (MRI) is often incorporated in a multimodal BLT system. In a previous study, the reconstruction of BLT based on the heterogeneous model with appropriate optical parameters provided more accurate results in localization and quantification than the homogeneous model [7]. Compared with MRI, the CT device has low cost and high resolution, thus many multimodal BLT systems use CT to provide the anatomic structure [12].
The heterogeneous BLT reconstruction usually involves four main steps. First, the anatomical structure of the imaged animal is obtained via main organ segmentation of the CT volume. Second, the segmented organs are discretized into finite elements such as tetrahedral or hexahedral meshes which are suitable for further reconstruction using finite element methods. Third, the two dimensional (2D) bioluminescence observation is mapped onto the 3D surface of the animal. Finally, the BLT problem is linearized into the form of matrix multiplication and solved using optimization methods. Organ segmentation plays an important role in the process of heterogeneous BLT reconstruction. However, the CT image is difficult to segment accurately because of the low contrast of soft tissues. The commonly used automatic segmentation methods are not appropriate for segmentation of low contrast structures in CT images, while manual segmentation is highly labor intensive, particularly with larger data sets, and prone to human error or biases. The use of intravenous contrast agents can improve the contrast of soft tissues [13], but it also increases the study cost and complexity. In order to reduce the processing complexity of segmentation, many studies focus on the estimation of the internal organs with atlas-based registration methods instead of accurate segmentation, which have gained success in human image segmentation [14][15][16]. Relatively fewer studies are concerned with the problem of small animal segmentation with atlas-based registration. Martin et al. presented a fully automated method for atlas-based whole-body segmentation in non-contrast enhanced Micro-CT data of mice [17]. But the single atlasbased registration is not accurate enough for estimating the internal organs because of the individual variances. Wang et al. constructed a statistical mouse atlas based on the statistical shape model and applied it to register with Micro-CT images for estimation of the mouse organ locations [18], whereas the mouse surface and organs of high contrast must be provided in their work. In order to relieve the data processing procedure while maintaining the registration accuracy, we propose a multi-atlas registration method to estimate the internal organs of a mouse in this study.
The reconstruction of BLT needs an accurate forward model of light transport in biological tissues. The Monte Carlo method [19-21] is regarded as the gold standard for simulation of light propagation in tissues. But it is time consuming because a large number of photons are required for reliable simulation results. The radiative transport equation (RTE) is usually used to model the forward problem of BLT [22]. However, RTE is difficult to solve numerically especially for complex anatomical structures. Hence the diffuse equations (DE) [23, 24] and simplified spherical harmonic equations (SPN) [25] are often used as the approximations for solving RTE. In this paper, we use DE approximation of RTE as the forward model for BLT. In order to solve the problem of light propagation in a complex structure with the finite element method (FEM) [26,27], the organs are often discretized into tetrahedral meshes. The generation of the tetrahedral meshes is a complex procedure, which can be implemented by using some commercial software such as Amira software (AMIRA, Mercury Computer Systems, Berlin, Germany) and Comsol Multiphysics (COMSOL Inc., Palo Alto, CA). Due to the expensive prices, this commercial software is not always available in most biological research laboratories. In this paper, we discretize the segmented volume of organs into adaptive hexahedral voxels. Compared with tetrahedral elements, hexahedral elements are easier to generate from CT volume and are more accurate in some research using the finite element analysis [28,29]. Moreover, the tetrahedral mesh has about five times more elements than the hexahedral one when they have the same number of nodes [30]. Therefore, the computational effort can be reduced if the volume is meshed using hexahedral elements.
The contribution of this paper can be summarized as follows. First, a multi-atlas registration method is conducted to estimate the internal organs of the imaged animal. The registration result shows that the multi-atlas registration is capable of accurately segmenting the organs with high efficiency. Second, an adaptive hexahedral voxel generation method is presented for solving the BLT problem with FEM. The merits are that it saves storage space and can be processed much easier than a tetrahedral mesh. These contributions alleviate the complexity of the data processing procedure of BLT while keeping an acceptable accuracy.

Forward and inverse problems of BLT using the hybrid optical/CT system
In this paper, a hybrid BLT system combined with optical/CT for small animal (a mouse in this paper) imaging is used for the study of the BLT problem. As depicted in Fig. 1, an X-ray tube and an X-ray detector are located on both sides of the animal holder. The rotational plate is perpendicular to the animal holder and rotates 360 degrees around it to capture the X-ray projection images. The anatomical structure of the mouse is obtained by using a cone-beam Feldkamp-Davis-Kress (FDK) reconstruction algorithm on the graphics processing unit using an acceleration scheme [31]. The bioluminescent image is captured by a scientific charge coupled device (CCD) camera whose axial direction is vertical to the X-ray central projection direction. The flow chart of BLT reconstruction is also shown Fig. 1. First, the top view bioluminescent image is captured by the CCD camera. Second, the 3D anatomical structure of the mouse is obtained by the CT device. The distributions of the organs are estimated via the multi-atlas registration method. Then, the segmented organs are discretized into adaptive hexahedral voxels and the 2D bioluminescent image is mapped onto the 3D surface of the mouse [32]. With the discretized voxels and light density on the boundary surface, the BLT problem can be linearized into the form of matrix multiplication by using the finite element method. Finally, the light source is reconstructed by optimizing the objective function of BLT based on Tikhonov regularization [33].
where Φ is the flux density, Q is the light source, Ω and ∂Ω represent the region of the object and its boundary respectively, a μ is the absorption coefficient, s μ is the scattering coefficient, and v is the unit outward normal on ∂Ω . , R x n n A x n n R x n n where n is the refractive index with Ω and n ′ is the refractive index in the surrounding medium. The surrounding medium is usually air. Therefore, ( , , ) R x n n ′ can be calculated by Since BLT is an ill-posed problem, the optimization objective function is defined by using Tikhonov regularization for the reconstruction: where sup S is the upper bounds of the source density, Λ stands for the weight matrix, λ is the regularization parameter, ( ) η ⋅ represents the penalty function, and meas Φ is the measured signal on the boundary. A modified Newton method and active set strategy is adopted to solve the minimum problem.

Adaptive voxel generation of the mouse
The adaptive hexahedral voxel is generated from the segmented CT data of the mouse. Hexahedral voxels of a fixed size cannot model the boundary of the organ accurately. So we generate hexahedral voxels of different sizes on the whole mouse body using a coarseto-fine strategy. As shown in Fig. 2, the largest hexahedral voxels termed as the first-level voxels are generated and labeled with the organ index according to its region. If a firstlevel voxel is on the boundary of different organs, the voxel is subdivided into 8 subvoxels and termed as the second-level voxel. Then, the second-level voxels are tested with the boundary condition and refined until each sub-voxel is completely within an individual organ or its size reaches the minimum size, i.e., the CT resolution. The boundary condition is defined with the organ consistency of the CT voxels overlapped by the sub-voxel. Suppose a sub-voxel overlaps n CT voxels. If all of the CT voxels are labeled with the same organ index, then the sub-voxel is defined to be located inside of an organ. Otherwise it is on the boundary of organs. The proposed adaptive generation method of hexahedral voxels is more effective in space partition and can reduce storage space compared with other types of meshes.  A cone-beamed Micro-CT device is utilized to obtain the CT volume data of the mouse. The voltage of the X-ray tube is set as 50 kilovolts and the acquisition time of each projection is set as 0.75 seconds for CT imaging. The reconstructed CT volume has a resolution of 0.409mm × 0.409mm × 0.409mm for each voxel and a dimension of 256 × 256 × 256. The mouse was anesthetized with isoflurane liquid and placed on the animal holder during the data acquisition. We constructed 50 mouse atlases from the acquired CT data. The organs of each mouse were segmented manually from the CT volume data. All of the segmented results were normalized by aligning all of the gravity centers to the position of (0, 0, 0) and rotating all of the volumes in the principal direction calculated by the principle component analysis (PCA). Some examples of the normalized mouse atlases are shown in Fig. 3.

Multi-atlas based mouse registration
Based on the above mouse atlases, the organ locations of a new CT image can be estimated by using the multi-atlas registration method. The posture and size of the mouse have a serious impact on registration accuracy. We define a similarity Sim to evaluate these factors as follows: where axial D is the axial distance between the observed mouse and atlas, which is defined as: where Si I is the center of the i th section along the principal direction of the observed mouse, and Ri I is that of the atlas. Dice is the Dice coefficient, which is calculated as: where R R and S R represent the organ regions of registration and segmentation respectively;  denotes the number of voxels. In Eq. 5) Cumulate the registration results and calculate the probable distribution of the segmented organs. Finally, the internal organs of the mouse can be segmented from the fusion of the registration results. Since each organ of the fusion volume is a connected region, we chose an adaptive seed-growing segmentation algorithm to label the organs.

Evaluation of the registration accuracy
The first experiment is conducted to study the accuracy of multi-atlas registration. A manually segmented mouse served as the test sample. Figure 4(a) and 4(b) show the surface and internal organs of the segmented mouse respectively. The surface of the mouse can be segmented easily by using the intensity thresholding method. Then, the multiple atlases are registered with the surface mask. Figure 4(c) shows the final distribution of the estimated organs after fusing the registration results. The segmented organs derived by using the adaptive seed-growing method are shown in Fig. 4(d). The single-based registration is also conducted for comparison. We also use the Digimouse [34] as the atlas to register with the test sample. Figure 4(e) shows the structure of the Digimouse for registration. The segmented organs based on the Digimouse registration are shown in Fig. 4(f). The mean Dice coefficients derived by the 10 selected atlases registration M Dice , the Dice coefficients M Dice ′ of the final result in Fig. 4(d), and the Dice coefficients S Dice of the single atlas registration result in Fig. 4(f) are compared in Table  1. The registration results demonstrate that the proposed multi-atlas registration method yields larger Dice coefficients than the single atlas based method, which implies that the proposed method can estimate the internal organs with high accuracy.

Numerical simulation of BLT
In this experiment, a digital mouse with an embedded light source is engaged in the simulation of BLT. The internal organs of the mouse, including bone, liver, heart, lungs, and kidneys are segmented manually from the Micro-CT volume data. The optical parameters of each organ under the wavelength of 620nm are listed in Table 2, including the absorption coefficient a μ , scattering coefficient s μ , anisotropy factor g, and refraction index n. The optical parameters are adopted from [35]. Monte Carlo simulation of BLT is conducted to get the forward result of BLT, which is used as the detected boundary condition for reconstruction. We modified the tetrahedron-based Monte Carlo simulator [20] to the hexahedron-based Monte Carlo simulator for simulation of BLT. The whole region of the mouse is discretized into the hexahedral voxels which consist of 6588 nodes and 4908 hexahedral elements. Figure 5

In vivo study of BLT
In this experiment, the fluorescent colloid containing the CT contrast agent was injected into the back of the mouse to model the bioluminescent source. The colloid can be segmented easily from the CT volume. In Fig. 6(a), the bioluminescent image is acquired by a highly sensitive CCD camera (Princeton Instruments PIXIS 2048B, Roper scientific, Trenton, NJ). A 40nm band-pass filter centered at 620nm is used to allow light transmission. The exposure time of the camera is set as 1.5 seconds for image acquisition. The organs of the CT volume are segmented by using the proposed multi-atlases registration method, which is shown in Fig. 6(b). The optical parameters of the organs under the 620nm wavelength are the same as that of the simulation study. In Fig. 6(c), the 2D bioluminescent image is mapped onto the mouse surface [32] which is discretized into hexahedral meshes. The volume rendering of the reconstructed result is given in Fig. 6(d).
It shows that the reconstructed light source colored in blue and the real one colored in red overlap together. The sectional views of the volume where the light source is located are presented in Fig. 6(e). The yellow and green dashed lines indicate that the reconstructed light source center and the real one are in good agreement. The result presents a 2.75mm distance between the position of the maximum reconstructed value and that of the real one, which suggests that the proposed method has good potential for practical applications.

Discussion
In this paper, the atlas registration algorithm is adopted to estimate the internal organ distributions of the mouse, which provides anatomical information to overcome the illposedness of BLT. The single atlas registration method is heavily affected by atlas selection. If the registered mouse's shape differs greatly from the selected atlas, it will cause great deviation and lead to a large registration error. Since the multiple atlases may involve mice with different features, multi-atlas based registration can handle the diversity of the observed mouse better. The experimental results in this paper have demonstrated the advantage of multi-atlas registration over single atlas registration. It should be noted that this study constructed mouse atlases only using healthy mice. Thus, the atlases are not suitable for registration of tumor bearing models, especially if the tumor is located deep inside the mouse's body. The problem can be tackled by detecting abnormal tissues first and then using this as a priori information for registration. This will be studied in the future.
The problem of BLT is linearized and solved by FEM based on the hexahedral voxel in this paper. The digital mouse model in section 3.2 has 5 times the tetrahedral elements when it is converted to tetrahedral meshes using the same number of nodes of the hexahedral voxels, which would cause a large computational burden. Therefore, BLT based on the hexahedral voxel is more efficient than the one based on the tetrahedral meshes. However, the comparison study of accuracy between BLT based on hexahedral voxels and tetrahedral meshes is not conducted in this paper. Extensive research will be performed in our future work.
The reconstruction of BLT is transformed to a least squares problem incorporated with Tikhonov regularization which uses the l 2 norms constraint. The result of BLT reconstruction indicated that the method can reconstruct the light source with acceptable precision. However, there are many other regularization techniques such as the l p norm regularization [36] and l 1 norm regularization [37]. Previous studies showed that these methods have effectiveness and potential in the application of quantitative BLT. We will test these regularization methods for BLT based on the hexahedral voxel in further research.

Conclusion
We present a fast reconstruction method for BLT. The aim of this method is to reduce the complexity of the data processing procedure for BLT by using a hybrid Optical/CT system. We adapt the multi-atlas registration algorithm to estimate the internal organs of the mouse instead of manual segmentation. The registration results indicate that the proposed method based on multi-atlas registration can estimate the distribution of organs with great efficiency and high accuracy. The hexahedral voxels are adaptively generated from the segmented CT volume. The adaptive hexahedral voxels for BLT reconstruction cannot only save the storage space but are also more efficient in data processing compared with tetrahedral meshes. The reconstructed results demonstrate that the proposed method can be employed to reconstruct the bioluminescent source effectively and accurately.