Open Access
Add to BibSonomy
Abstract
To alleviate the ill-posedness of bioluminescence tomography (BLT) reconstruction, anatomical information from computed tomography (CT) or magnetic resonance imaging (MRI) is usually adopted to improve the reconstruction quality. With the anatomical information, different organs could be segmented and assigned with appropriate optical parameters, and the reconstruction could be confined into certain organs. However, image segmentation is a time-consuming and challenging work, especially for the low-contrast organs. In this paper, we present a BLT reconstruction method in conjunction with an organ probability map to effectively incorporate the anatomical information. Instead of using a segmentation with a fixed organ map, an organ probability map is established by registering the CT image of the mouse to the statistical mouse atlas with the constraints of the mouse surface and high-contrast organs (bone and lung). Then the organ probability map of the low-contrast organs, such as the liver and kidney, is determined automatically. After discretization of the mouse torso, a heterogeneous model is established as the input for reconstruction, in which the optical parameter of each node is calculated according to the organ probability map. To take the advantage of the sparse Bayesian Learning (SBL) method in recovering block sparse signals in inverse problems, which is common in BLT applications where the target distribution has the characteristic of sparsity and block structure, a two-step method in conjunction with the organ probability map is presented. In the first step, a fast sparse algorithm, L1-LS, is used to reveal the source distribution on the organ level. In the second step, the bioluminescent source is reconstructed on the pixel level based on the SBL method. Both simulation and in vivo experiments are conducted, and the results demonstrate that the organ probability map in conjunction with the proposed two-step BLT reconstruction method is feasible to accurately reconstruct the localization of the bioluminescent light source.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
Corrections
Wanzhou Yin, Xiang Li, Qian Cao, Hongkai Wang, and Bin Zhang, "Bioluminescence tomography reconstruction in conjunction with an organ probability map as an anatomical reference: erratum," Biomed. Opt. Express 14, 1243-1244 (2023)https://opg.optica.org/boe/abstract.cfm?uri=boe-14-3-1243
1. Introduction
Bioluminescence tomography (BLT) is a rapidly developing optical molecular imaging modality, which can be used to monitor physiological processes on molecular and cellular levels [1,2]. BLT has the advantages of no ionizing radiation, in vivo and high sensitivity, and provides a powerful tool for tumor monitoring, drug development and disease occurrence mechanism [3–5].
In bioluminescence imaging (BLI), the 2D bioluminescence signal on the tissue surface is captured directly. However, due to the strong scattering and absorption of biological tissue, the position of the internal light source cannot be correctly determined by the measurements on the surface. In contrast, the BLT could show the bioluminescent targets inside the tissue in 3D using the light signals (either single or multiple spectral data) measured on the tissue surface by adopting a sophisticated image reconstruction algorithm. By comparing the predicted light intensity on the tissue surface with the measured intensity and iteratively updating the unknown bioluminescent source distribution, the reconstruction algorithm finally reveals the target position [2,5,6].
The predicted light intensity is determined with a light propagation model which requires optical properties and the geometry of the animal. The animal’s geometry for modeling light propagation can be obtained either by using structured light for capturing the animal’s surface or by other imaging modalities, such as x-ray computed tomography (CT) or magnetic resonance imaging (MRI) [6–11]. The former approach captures a single view of the animal surface which leads to a homogeneous model, while the later approach utilizes a CT or MRI volume which is segmented to organs. Then a heterogeneous model could be built based on the segmented high-resolution CT or MRI anatomical images with each organ is assigned with appropriated optical property.
Although the heterogeneous model is more accurate than the homogeneous one in modeling light propagation for BLT, the organ segmentation is time-consuming and challenging, especially for the low-contrast CT images. To alleviate the burden of organ segmentation from CT and facilitate BLT reconstruction, registration based on atlas has been adopted and attained researcher’s attention [12–15]. Ren et al. [16] presented a registration method to estimate the organ’s location of the imaged mouse based on multi-atlas. The atlas is constructed using CT datasets from 50 mice which are manually segmented into different organs respectively. Compared with single atlas-based registration, the multi-atlas-based registration can reduce individual variance and get more accurate registration. Recently, Klose et al. [17] proposed a computer-aided analysis method for BLT reconstruction. They innovated a body-conforming animal mold to spatially scale the animal, and then using the mouse atlas to register the animal constrained in the mold. In combination with a pre-solved multispectral light propagation model of each mold, the in vivo biodistribution of bioluminescent reporters can be rapidly calculated. The results demonstrated that the BLT in combination with a statistical mouse atlas could facilitate quantitative analysis of bioluminescence images.
Different from the works of [17], our group have proposed a more general statistical mouse atlas for organ segmentation which does not need an animal mold. The statistical mouse atlas was generated from 103 mice [18,19]. In practical, the CT image of the mouse is first registered to the statistical mouse atlas with the constraints of mouse surface and high-contrast organs (such as the bone and lung). Then the low-contrast organs, such as the liver and kidney, are estimated automatically through the statistical mouse atlas model. In our recent work [20], with the estimated shapes of the organs from the statistical mouse atlas, we have demonstrated the effectiveness in promoting BLT reconstruction for the multimodality BLT/CT system, and the results show that the bioluminescent sources imbedded in the liver and lung can be accurately reconstructed based on the statistical mouse atlas, which is comparative to the reconstruction result using an accurate organ segmentation.
In this paper we presented a two-step BLT reconstruction method in conjunction with an organ probability map to incorporate the anatomical information. Different from our previous work [20], where the shapes of the low-contrast organs are fixed and calculated based on the mean shape of the statistical mouse atlas, in this paper, the low-contrast organs are estimated as an organ probability map. The optical property map of the mouse is established by weighting the optical properties of different organs based on the organ probability map, and it is served as an input for BLT reconstruction. To take the advantage of sparsity in BLT reconstruction, the first-step reveals the distribution of the bioluminescent source on the organ level by utilizing a fast sparse algorithm, L1-LS [21,22]. The sparse Bayesian Learning (SBL) method [23–26] has the advantage of recovering block sparse signals in inverse problems, which is suitable for the applications in BLT where the target distribution has the characteristic of sparsity and block structure. The SBL could get better performance if the block structure knowledge is known in advance. So, a two-step reconstruction method is proposed, and in the second-step, the bioluminescent source is reconstructed on the pixel level based on the SBL method with the initial result from L1-LS. Simulation and in vivo experiments are conducted to evaluate the feasibility of the proposed method to localize the bioluminescent sources in different organs.
2. Materials and method
2.1 BLT reconstruction
The radiative transfer equation (RTE) has been proved as the most comprehensive model for describing light transportation in biological tissue. But it is challenging to find the numerical solution of RTE, especially for the object with irregular shape. So many approximate methods to the RTE have been developed to overcome the difficulties of solving the RTE. The diffusion equation (DE) which is the 1st order approximation of the RTE, is widely employed to model light transportation in biological tissue. However, it is only valid under the condition that the reduced scattering coefficient is much larger than the absorption coefficient (${\mu _s}^{\prime}/({\mu _s}^{\prime} + {\mu _a}) \approx 1$). For wavelength shorter than 620nm, the solution of DE is less accurate because the light absorption becomes stronger (${\mu _s}^{\prime}/({\mu _s}^{\prime} + {\mu _a}) < 1$) [6]. In this paper, the 3rd order simplified spherical harmonics function (SP3) to the RTE is used to model light propagation in biological tissues [6,27,28]. As a high order approximation to the RTE, SP3 can balance the contradiction between accuracy and computational effort comparing with DE. The steady state SP3 equation is expressed as follow:
If the source is a point source, then its 3D photon flux distribution can be expressed as Green’s function ${G^{em}}({r_d},r)\textrm{ = }{\varphi _1}(r) - \frac{2}{3}{\varphi _2}(r)$ using the SP3 method. According to the 1st order Born approximation, the measurement on the surface point ${r_d}$ can be expressed as
where $\Theta $ is the gain factor related to the system settings, and $x(r)$ is the unknown bioluminescent source distribution.Assume the imaged object is discretized into N mesh nodes, and then the optical parameter of each node is assigned through a heterogeneous optical property map. With an appropriate boundary condition, the Green’s function of Eq. (1) can be solved by finite element method (FEM). The linear relationship between the measurements on the object surface, y, and the unknown bioluminescent source distribution, x, can be established as
where W represents the system matrix with the size of ${N_d} \times {N_v}$, ${N_d}$ is the number of measurements on the surface, and ${N_v}$ is the number of the discretized nodes. Here, multispectral measurements are acquired to improve BLT reconstruction, and the relative spectral weights caused by the source spectrum, filters and the CCD quantum efficiency at different wavelengths are also considered into Eq. (3). A modified version of the open source NIRFAST software is used to calculate W [29,30]. Due to the measurements are much fewer than the unknowns (${N_d} \ll {N_v}$), the BLT reconstruction is ill-posed and underdetermined, so it is impossible to solve Eq. (3) directly.2.2 Statistical shape model
Our group had previously constructed a statistical mouse atlas based on 103 manually segmented mouse CT images [18–20]. In actual operation, the mouse is first scanned for CT imaging followed by an automated organ segmentation. The low-contrast organs, such as the liver and kidneys, can be automatically estimated by registering the high-contrast organs (mouse surface, lung and bone) to the statistical mouse atlas. Then the estimated organs providing structure information is used as input for BLT reconstruction.
Different from our previous work [20] where the shapes of the low-contrast organs are estimated based on the mean value of each organ of the statistical mouse atlas, in this paper, the low-contrast organs are established as an organ probability map, as shown in Fig. 1. In the organ probability map, the position far away from the center of the organ has lower probability. The organ probability map of the mouse is more appropriate for compensating anatomical variations comparing with the fixed mean shape.
Fig. 1. (a) The organ probability map estimated from the statistical shape model. The probability distribution of (b) the liver, (c) the right kidney and (d) the soft tissue.
The organ probability map is then translated into an optical property map, including the absorption coefficient ${\mu _a}$ and the scattering coefficient ${\mu _s}$, which contain the different optical properties of organs [31,32], as shown in Table 1. The optical property maps of ${\mu _a}(\lambda )$ and $\mu _s^{\prime}(\lambda )$ are built for four wavelengths between 590nm and 650nm with 20 nm interval. With the relationship ${\mu ^{\prime}_s}\textrm{ = }{\mu _s}(1 - g)$, the isotropic parameter g is set to 0.9 in all the simulation and in vivo experiments, then the scattering coefficient ${\mu _s}$ could be easily calculated.
Table 1. Optical parameters of different regions under four wavelengths [31] ($\textrm{m}{\textrm{m}^{ - 1}}$)
The weighted optical parameter of each mesh node is calculated as
where $\mu _i^j(\lambda )$ is the optical parameter of the jth organ (j = 1, 2, …, Nseg), representing absorption coefficient ${\mu _a}(\lambda )$ or scattering coefficient ${\mu _s}(\lambda )$, $p_i^j$ is the organ probability of organ j at position I, and $\sum\limits_j^{Nseg} {p_i^j} = 1$. The optical property maps were then input to the light propagation model of the BLT reconstruction.2.3 Two-step BLT reconstruction
Taking advantage of the organ probability map, here we developed a two-step reconstruction method for BLT reconstruction. The first step is to calculate the distribution of the bioluminescent source at organ level. In the second step, a high-resolution solution at voxel level is reconstructed based on the prior information from the first step using the sparse Bayesian learning method. The details of the two-step reconstruction method are introduced as follows.
Step 1: BLT reconstruction on organ level. Through registration with the statistical mouse atlas, as described in Section 2.2, the imaged mouse is segmented into Nseg organs. Then a segmentation matrix G (${N_\textrm{v}}\mathrm{\ \times }{N_{\textrm{seg}}}$) is constructed, where the element ${G_{ij}}$ is the probability of jth organ at mesh grid i. Assuming the signal intensities are equivalent inside a specific organ, then a low dimensional linear equation can be subsequently constructed:
where $\bar{W} = W \cdot G$ is the new system matrix with the size of ${N_\textrm{d}} \times {N_{\textrm{seg}}}$, the unknown vector $\bar{x}$ includes Nseg elements representing the average intensity of the Nseg organs. Equation (5) can be converted into an optimization problem where $\lambda$>0 is the regularization parameter. Here, Eq. (6) is solved using the widespread sparse based lease square L1-LS algorithm with the constraint of non-negativity. The L1-LS algorithm is developed based on an interior-point method for L1-resularized least squares problems and is fast for large sparse problems. Owning to the small dimension of $\bar{W}$, the solution could be very fast. Consequently, the solutions for all mesh nodes are obtained through $x = G\bar{x}$, and only the positions with non-zero solutions are selected for the 2ed step.Step 2: High-resolution BLT reconstruction. In most applications, the bioluminescence source localizes only a small portion of the whole field of view, i.e., the distribution has the characteristic of sparsity and block structure. Aiming at this feature, Zhang et al. [33] proposed a sparse Bayesian learning (SBL) algorithm for recovering block sparse signals. The SBL algorithm could take advantage of the initial results at organ level reconstructed from the first step. Here, we adopted the SBL algorithm on the high-resolution BLT reconstruction.
For the sake of explanation, a brief introduction of SBL is listed as follows [23–25,33]. Assume the unknowns ${x_i}$ in the i-th organ satisfying a multivariate Gaussian distribution:
where ${\gamma _i}$ is a non-negative parameter that control the sparsity of ${x_i}$,when ${\gamma _i}\textrm{ = }0$, the corresponding ${x_i}$ become zero. ${B_i}$ is a positive definite matrix, representing the correlation between elements in ${x_i}$. Assuming that the unknowns ${x_i}$ are mutually independent, then the prior of unknown x in Nk organs is: where ${C_0} = diag({\gamma _1}{B_1},{\gamma _2}{B_2}, \cdots ,{\gamma _g}{B_g})$. ${B_i}$ is an all-one matrix if all elements in ${x_i}$ take the same value.Under the assumption that the noise v is a white noise and satisfies $p(v;\lambda )\sim N(0,\lambda I)$, where $\lambda$ is a positive scalar, the conditional probability is given by
Based on the Bayesian rule, the posterior of x is obtained through
withWhen parameters $\lambda ,\gamma$ and $B$ are estimated, ${w_x}$ can be regarded as an estimation of x.
A type II maximum likelihood is used to estimate the parameters $\lambda ,\gamma$ and B, which leads to minimizing the following cost function:
Taking the derivative of Q with respect to $\lambda ,\gamma$ and B respectively, the learning rules can be derived as follows:
Substitute the estimated parameters into Eq. (10), and we can obtain the estimation of the solution x. The flow chart of the algorithm is shown in Fig. 2.
2.4 Evaluation parameters
To quantitatively evaluate the performance of the proposed BLT reconstruction method, the following three parameters have been calculated.
- (1) Center of mass (COM) and 3D offset (Doffset) between the reconstructed source and the true source. where (x, y, z) and (x0, y0, z0) are the COM’s coordinates of the reconstructed and the true sources, respectively. The Doffset represents the position accuracy of reconstruction.
- (2) In order to analyze the shape similarity of the reconstructed and the true sources, the Dice coefficient is calculated. The Dice coefficient is given by: where A and B are the volumes of the reconstructed and true source, respectively. After reconstruction, the voxels with intensity higher than 10% of the maximum of the intensity are regarded as the reconstructed volume A. With the known true source volume B in simulations, the Dice coefficient could be calculated.
- (3) Contrast to noise ratio (CNR) is used to evaluated the quality of reconstructed results. Given a region of interest (ROI), the CNR is defined as follow:$$CNR = \frac{{{\mu _{ROI}} - {\mu _{BCK}}}}{{{{({\omega _{ROI}}\sigma _{ROI}^2 - {\omega _{BCK}}\sigma _{BCK}^2)}^{1/2}}}}, $$where ${\mu _{ROI}}$ and ${\mu _{BCK}}$ are the mean intensity values, $\sigma _{ROI}^2$ and $\sigma _{BCK}^2$ are the variances, and ${\omega _{ROI}}$ and ${\omega _{BCK}}$ are the weighing factors of the ROI volume and background volume, respectively. After reconstruction, the reconstructed value is converted to the range of [0, 1] through normalizing to the maximum value. Then the voxels with intensity higher than the threshold of 0.1 is regarded as the reconstructed BLT volume. After subtracting the reconstructed BLT volume, the rest of the volume is regarded as the background volume. The higher value of CNR indicates better quality of the reconstruction
3. Experiments and results
In order to assess the feasibility of the proposed two-step BLT reconstruction method in conjunction with the statistical shape model, simulations and in vivo experiments were conducted, and the results were compared with the classical sparse reconstruction method, L1-LS. Additionally, the reconstruction performances were also compared with different organ segmentation methods, including the accurate organ segmentation (manual segmentation), the fixed mean shape segmentation, and the proposed organ probability map.
3.1 Simulation configuration
The 3D mouse model, Digimouse [34], is used for simulation. Performance of the proposed reconstruction method is test by artificially placing tumor targets in different organs. The heart, lung, liver, kidney, bone and soft tissue of the Digimouse are selected and regarded as the true organ localization, as shown in Fig. 3(a). The organ probability map is established through registration of the Digimouse to the statistical shape model. The mouse surface, bone and lung are registered by fitting the statistical shape model for their high-contrast in CT image, then the organ probability maps, such as the heart, liver and kidney are subsequently estimated from the statistical shape model, as illustrated in Figs. 3(b-d).
Fig. 3. (a) The heterogeneous model of the Digimouse. (b∼d) Registration of the organ probability map to the Digimouse. (e) The discretized mesh for generating multispectral measurements.
To save computational efforts, only the torso part between the dish line, as shown in Fig. 3(a), is cropped for meshing. Multispectral data at 590, 610, 630 and 650 nm are acquired for the purpose of increasing the surface measurements and promoting the BLT reconstruction. To make the simulation close to the real situation, the multispectral measurements are added with 5% of white noises.
To test the performance of the proposed algorithm on organs with different optical parameters, four simulation experiments are conducted with the sources located in the liver, kidney, abdomen and lung, respectively. The simulated multispectral measurements are calculated by the SP3 method mentioned in Eq. (1) on a fine mesh discretized from the mouse torso with true organ localization. As shown in Fig. 3(e), the mesh has ∼290,000 tetrahedron elements and ∼52,000 nodes, and the mean element edge size is 0.6mm. Then, the simulated multispectral measurements and the cropped mouse torso with the organ probability map are used in the following reconstruction.
3.2 Simulation 1: light source in the liver
In this simulation, a sphere light source with a radius of 1.5 mm and the center coordinate of (16.6, 49.8, 8.3) mm is placed in the liver as shown in Figs. 4(a1∼a3). The mouse torso is discretized into a mesh (∼150,000 tetrahedron elements and ∼27,000 nodes, and the mean element edge size is 0.8mm) and the optical parameters of each node are assigned based on the organ probability map as shown in Section 2.2. The mesh, the optical parameters and the simulated multispectral measurements (calculated in Section 3.1) are used for BLT reconstruction.
Fig. 4. Reconstruction results of a single source in the liver. (a1) 3D rendering of the true source distribution, and (a2, a3) the transverse and coronal views of the true source overlapped with the CT images. (b1) 3D rendering of the L1-LS reconstruction result and (b2, b3) the reconstruction result overlapped with the CT images. (c1) 3D rendering of the two-step reconstruction result with the fixed mean shape segmentation, and (c2, c3) the reconstruction result overlapped with the CT images. (d1) 3D rendering of the two-step reconstruction result with the proposed organ probability map, and (d2, d3) the reconstruction result overlapped with the CT images. All the BLT results are normalized to the maximum value.
After reconstruction using the first-step of the proposed method, the normalized intensities of each organ are: heart 0.014, liver 0.071, right kidney 0, left kidney 1.0, lung 0.026, born 0 and soft tissue 0. So only the non-zero intensity organs (the heart, liver, left kidney and lung) are selected for the second-step high-resolution reconstruction, and the final reconstructed results are shown in Fig. 4(d1). The evaluation parameters have been shown in Table 2, the COM coordinate is (16.8, 50.0, 8.0) mm, Doffset is 0.4mm, Dice is 0.71 and CNR is 81.7. The fusion of BLT reconstruction and CT image is shown in Figs. 4(d2-d3). In contrast, the BLT reconstruction with the same configuration but the L1-LS method is also performed, and the corresponding results are shown in Figs. 4(b1-b3). The evaluation parameters of the L1-LS method are also listed in Table 2 for comparing.
Table 2. Quantitative comparison of the reconstruction results based on the proposed two-step method and the L1-LS method
The BLT reconstruction with the two-step method but using the fixed mean shape segmentation is shown in Figs.4(c1-c3). The quantitative comparison of reconstruction results with different segmentation methods is listed in Table 3.
Table 3. Quantitative comparison of the reconstruction results with different segmentation methods
3.3 Simulation 2: light source in the kidney
In this simulation, a sphere light source (1.5 mm radius and center at (25.0, 41.1, 14.0) mm) is placed in the kidney, as shown in Figs. 5(a1-a3). The torso part is discretized into a mesh (∼170,000 tetrahedron elements and ∼31,000 nodes, and the mean element edge size is 0.8mm) and the optical parameters of each node are assigned based on the organ probability map as shown in Section 2.2. The mesh, the optical parameters and the simulated multispectral measurements are used for BLT reconstruction.
Fig. 5. Reconstruction results of a single source in the kidney. (a1) 3D rendering of the true source distribution, and (a2, a3) the transverse and coronal views of the true source overlapped with the CT images. (b1) 3D rendering of the L1-LS reconstruction result and (b2, b3) the reconstruction result overlapped with the CT images. (c1) 3D rendering of the two-step reconstruction result with the fixed mean shape segmentation, and (c2, c3) the reconstruction result overlapped with the CT images. (d1) 3D rendering of the two-step reconstruction result with the proposed organ probability map, and (d2, d3) the reconstruction result overlapped with the CT images. All the BLT results are normalized to the maximum value.
After reconstruction using the first-step of the proposed method, the normalized intensities of each organ are: heart 0, liver 0, right kidney 1.0, left kidney 0.267, lung 0, bone 0.006 and soft tissue 0. So only the non-zero intensity organs (the right kidney, left kidney and bone) are selected for the second-step high-resolution reconstruction, and the final reconstructed results are shown in Fig. 5(d1). The COM coordinate is (24.4, 41.3, 13.9) mm, Doffset is 0.7mm, and the Dice is 0.57 and CNR is 337.3. The fusion of BLT reconstruction and CT image is shown in Figs. 5(d2-d3).
In contrast, the BLT reconstruction with the same configuration but the L1-LS method is also performed, and the corresponding results are shown in Figs. 5(b1-b3). The evaluation parameters of L1-LS method are also listed in Table 2 for comparison.
The BLT reconstruction with the two-step method but using the fixed mean shape segmentation is shown in Figs.5(c1-c3). The quantitative comparison of reconstruction results with different segmentation methods is listed in Table 3.
3.4 Simulation 3: light source in the abdomen
In this simulation, a sphere light source is placed in the abdomen, as shown in Figs. 6(a1-a3). The torso part was discretized into a mesh (∼170,000 tetrahedron elements and ∼31,000 nodes, and the mean element edge size is 0.8mm) for BLT reconstruction.
Fig. 6. Reconstruction results of a single source in the abdomen. (a1) 3D rendering of the true source distribution, and (a2, a3) the transverse and coronal views of the true source overlapped with the CT images. (b1) 3D rendering of the L1-LS reconstruction result and (b2, b3) the reconstruction result overlapped with the CT images. (c1) 3D rendering of the two-step reconstruction result with the fixed mean shape segmentation, and (c2, c3) the reconstruction result overlapped with the CT images. (d1) 3D rendering of the two-step reconstruction result with the proposed organ probability map, and (d2, d3) the reconstruction result overlapped with the CT images. All the BLT results are normalized to the maximum value.
After reconstruction using the first-step of the proposed method, the normalized intensities of each organ are: heart 0, liver 0, right kidney 1.0, left kidney 0, lung 0, bone 0 and soft tissue 0.13. Then only the non-zero intensity organs (the right kidney and abdomen) are selected for the second-step high-resolution reconstruction, and the final reconstructed results are shown in Fig. 6(d1). The COM coordinate is (9.6, 32.7, 7.8) mm, Doffset is 0.3mm, the Dice is 0.72 and the CNR is 103.7. The fusion of BLT reconstruction and CT image is shown in Figs. 6(d2-d3). Figure 6(b1) is the BLT reconstruction based on the L1-LS method, and the corresponding fusion of BLT reconstruction and CT image is shown in Figs. 6(b2-b3). The evaluation parameters of L1-LS method are also saved in Table 2 for comparison.
The BLT reconstruction with the two-step method but using the fixed mean shape segmentation is shown in Figs.6(c1-c3). The quantitative comparison of reconstruction results with different segmentation methods is listed in Table 3.
3.5 Simulation 4: double light sources in the lungs
In this simulation, two sphere light sources (1.5 mm radius and the centers are (16.7, 56.9, 13.0) mm and (23.2, 56.9, 11.0) mm) are placed in the lungs, as shown in Figs. 7(a1-a5). The torso part was discretized into a mesh of ∼110,000 tetrahedron elements and ∼21,000 nodes for BLT reconstruction with the mean element edge size is 0.8mm. With the proposed reconstruction method, the mean intensities of each organ after the first-step are: heart 0, liver 1.74e-5, right kidney 0, left kidney 1.0, lung 8.74e-5, bone 0 and soft tissue 1.77e-6. Then only the non-zero intensity organs (the liver, left kidney, lungs and the soft tissue) are selected for the second-step high-resolution reconstruction, and the final reconstructed results are shown in Fig. 7(d1). The fusion of BLT reconstruction and CT image is shown in Figs. 7(d2-d5). For comparison, the BLT reconstruction with the same configuration but the L1-LS method is also performed, and the corresponding results are shown in Figs. 7(b1-b5). The evaluation parameters of two methods were also save in Table 2 for comparison.
Fig. 7. Reconstruction results of double sources in the lungs. (a1) 3D rendering of the true source distribution, and (a2, a3) the transverse and coronal views of the true source overlapped with the CT images. (b1) 3D rendering of the L1-LS reconstruction result and (b2, b3) the reconstruction result overlapped with the CT images. (c1) 3D rendering of the two-step reconstruction result with the fixed mean shape segmentation, and (c2, c3) the reconstruction result overlapped with the CT images. (d1) 3D rendering of the two-step reconstruction result with the proposed organ probability map, and (d2, d3) the reconstruction result overlapped with the CT images. All the BLT results are normalized to the maximum value.
The BLT reconstruction with the two-step method but using the fixed mean shape segmentation is shown in Figs.7(c1-c5). The quantitative comparison of reconstruction results with different segmentation methods is listed in Table 3.
3.6 In vivo experiments
To further test the performance of the proposed method, a group of in vivo experiments have been conducted in the BLT/CBCT system [35]. A bioluminescence source (a cylinder with 0.9 mm diameter and 2 mm length, Trigalight, mb-microtecag, Switzerland) is placed in abdomen though surgery when the mouse is anesthetized. The mouse is taped on the animal stage in supine position for imaging. Multispectral measurements at 590, 610, 630, 650 nm are acquired with CCD camera. The CT data is collected by CBCT system, and the exact position of the light source can be localized because its high contrast in CT image.
Due to the low-contrast of soft tissue, the organs such as liver and kidney are difficult to distinguish. The mouse surface and the high contrast organs (the bone, heart and lung) are quickly segmented through thresholding. After registering with the statistical mouse atlas, the organ probability map is established using the method described in Section 2.2 and is shown in Figs. 8(a1-a3). After segmentation, the optical parameter of each organ is calculated with Eq. (4).
Fig. 8. Reconstruction results of the in vivo experiments. (a1-a3) Registration of the organ probability map to the mouse CT image. (b1) Surface optical signals collected during the experiment, and (b2) the surface optical signal distribution on the discretized mesh. (c1) 3D rendering of the reconstruction result using the proposed two-step method, and (c2, c3) the coronal and transverse views of the overlapped BLT and CT images. (d1) 3D rendering of the reconstruction result using the L1-LS method, and (d2, d3) the coronal and transverse views of the overlapped images. All the BLT results are normalized to the maximum value.
A section of the mouse torso is cropped from the CBCT image and is discretized into a 3D tetrahedral mesh (including ∼97,000 tetrahedron elements and ∼18,000 nodes, and the mean element edge size is 0.8mm) for BLT reconstruction. The multispectral BLIs are then mapped on the mesh surface, as shown in Fig. 8(b2). With the proposed two-step reconstruction method, the normalized mean intensities of each organ after the first-step are: heart 0, liver 0, right kidney 1.0, left kidney 0, lung 0, bone 0, abdomen 0.012. Then only the non-zero intensity organs (the right kidney and the soft tissue) are selected for the second-step high-resolution reconstruction, and the final reconstructed results are shown in Fig. 8(c1). The center of the source is (8.7, 28.2, 2.1) mm, and the COM of reconstructed source is (9.1, 27.9, 1.5) mm, so the Doffset is 0.84 mm. The fusion of BLT reconstruction and CT image is shown in Figs. 8(c2-c3). For comparison, the 3D distribution of the reconstructed source based on the L1-LS method is also calculated and shown in Fig. 8(d1). The corresponding COM is (9.2, 28.4, 3.1) mm, and the Doffset is 1.16 mm. The fusion of BLT reconstruction and CT image is shown in Figs. 8(d2-d3).
4. Discussion
BLT reconstruction in conjunction with heterogeneous structure information can provide more accurate results than the homogeneous model [6,36,37]. Micro-CT is widely adopted in the multimodal BLT system to provide anatomic structure, but the segmentation of soft tissue is a challenging problem for its low contrast in CT image when contrast enhancement is absent. In this paper, the organ probability map is established to facilitate soft tissue organ segmentation of the mouse CT images. By registering the statistical shape model to the high-contrast organs (surface, bone and lung) of the CT images, the low-contract organs could be effectively segmented, which can be served as an anatomical information for BLT reconstruction and alleviate the ill-posedness of inverse solution. The feasibility of the organ probability maps in conjunction with the BLT reconstruction method is estimated through simulation and in vivo experiments.
Different from our previous work [20] in which the shapes of the low-contrast organs are estimated using the mean value from the statistical mouse atlas, in this paper, the organ probability map of the low-contrast organ is established. Because our statistical mouse atlas is constructed using 103 mice of different sexes, strains, weights and postures, and in combination with the conditional Gaussian model, so it is charactered with the capability of capturing the inter-subject and inter-organ shape variations [18,19]. So, the organ probability map is more appropriate for compensating anatomical variations comparing with the fixed mean shape.
With the organ probability map, the optical property map of different organs is built as shown in Section 2.2 and Fig. 1. Therefore, the structural information in mouse CT image is efficiently adopted. Besides, we developed a two-step method for BLT reconstruction. In the first step, the structural information is utilized to confine the solution in a low dimensional domain and the mean intensity of each organ is reconstructed based on a classical sparse method L1-LS. In the second step, the high-resolution 3D distribution of the bioluminescence source is reconstructed using the sparse Bayesian learning method.
Four simulation experiments, from different position and number of the light source, were conducted to test the performance of the proposed method. As illustrated in Section 3.2-3.5, we can see that: (1) the mean Doffset of the proposed method is 0.32mm and it is significantly smaller than that of the L1-LS method (0.72 mm), which means that the 3D of the localization of the reconstructed source is more accurate based on the proposed two-step method than that of the classical L1-LS method; (2) the mean Dice of proposed two-step reconstruction method (0.70) is much larger than that of the L1-LS method (0.29), which means that the size of reconstructed source is closer to the true source based on the proposed two-step method than that of the classical L1-LS method; and (3) the mean CNR of reconstructed image based on the proposed two-step method (181.1) is much better than that of the L1-LS method (13.3).
Additionally, the BLT reconstruction performances with different organ segmentation methods, including the accurate organ segmentation (manual segmentation), the fixed mean shape segmentation and the proposed organ probability map, were shown in Table 3 and Figs. 4–7. From the quantitative comparison shown in Table 3, it is notable that the CoM and Dice of the reconstruction using the proposed organ probability map are closer to those using the accurate organ segmentation, and superior to those results reconstructed using the fixed mean shape.
Meanwhile, as shown in Fig. 4(b1), Fig. 5(b1), Fig. 6(b1) and Fig. 7(b1), the reconstructed images are contaminated with noticeable noise when reconstructed using the L1-LS method. In contrast, the proposed two-step reconstruction method can greatly suppress noise signal, as shown in Fig. 4(d1), Fig. 5(d1), Fig. 6(d1) and Fig. 7(d1). Additionally, in vivo experiments were carried out to verify the feasibility of the proposed method. As illustrated in Fig. 8, the performance of the proposed two-step reconstruction method is superior to the L1-LS method in terms of the size of reconstructed source, the 3D COM offset and the image quality, which is consistent with the observation of the simulation results. The simulation and in vivo experiments indicate that the proposed method can reconstruct the bioluminescent source more accurate than the classical L1-LS method. In simulations the kidney is always selected after the first-step reconstruction, even the source was far away from kidneys. This is mainly because the segmented kidney from the organ probability map has noticeable deviation from the ground truth. The mismatch of optical parameters would bias the first-step reconstruction. With the second-step, the target could be reconstructed accurately. The computation time for the simulation and mouse experiments was all less than 3 minutes for the two-step BLT reconstruction using a 64-bit desktop with an Intel i7-9700K CPU 3.60GHz processor and 32 GB of memory.
There are some aspects which are not considered in this paper. First, the shape of the source is concentrated as a small sphere which is a special case in real situation. Therefore, more experiments containing different shapes should be carried out to fully verify the effectiveness of the proposed algorithm in the future work.
Second, many other sparse regularization methods [38–41], such as L2,1 and Lp norm regularization, have shown great potential in effectively solving the ill-posed BLT reconstruction problem. The purpose of this paper is to show the effectiveness of organ probability map in conjunction with the two-step BLT reconstruction method, rather than comprehensive comparison, so only the classical L1-LS method is adopted for comparison. Additionally, it is possible to replace the algorithms used in the paper to achieve two-step reconstruction, especially the L1-LS used in the first step.
Besides, the deep learning related methods have been widely used and become dominant in field of medical imaging. Recently, researchers have utilized the multilayer perceptron [42] and one-dimensional convolutional neural network [43] to solve the BLT reconstruction problem, and the results are promising over the traditional methods. Taking the advantage of the deep learning strategy to further improve the BLT reconstruction quality will be concerned in the future.
Last, in the in vivo study, a small Trigalight light source is placed into the abdomen of the mouse to evaluate the performance of the proposed method, which has high contrast in CT image and can be located accurately as reference for algorithm development. We only analyzed the localization accuracy of the target, the quantification of the bioluminescence, which is essential in biological research as well, is not considered here and will be studied in our future work.
5. Conclusion
In this study, the organ probability map in conjunction with a two-step BLT reconstruction method is presented and estimated through simulation and in vivo experiments. The organ probability map is established by registering with the statistical mouse atlas, and could effectively reveal the low-contrast organs in CT images which are difficult to segment. The heterogeneous model with appropriate optical parameters is built for reconstruction. Then a two-step method is presented for reconstructing the 3D distribution of the bioluminescence source. The simulation and in vivo experiments have demonstrated the feasibility of the proposed algorithm. The developed two-step method in conjunction with the organ probability map for BLT reconstruction combines L1 regularization and SBL algorithm and takes the advantages of priori information, significantly alleviating the ill-posedness of BLT reconstruction and improving the quality of the results.
Funding
Natural Science Foundation of Liaoning Province (2019-ZD-0645); National Natural Science Foundation of China (61571076, 81971693); Fundamental Research Funds for the Central Universities (DUT19JC01, DUT20YG122); Dalian Engineering Research Center for Artificial Intelligence in Medical Imaging.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008). [CrossRef]
2. C. Qin, J. Feng, S. Zhu, X. Ma, J. Zhong, P. Wu, Z. Jin, and J. Tian, “Recent advances in bioluminescence tomography: methodology and system as well as application,” Laser Photonics Rev. 9, 94–114 (2015). [CrossRef]
3. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef]
4. M. Keyaerts, V. Caveliers, and T. Lahoutte, “Bioluminescence imaging: looking beyond the light,” Trends Mol Med 18(3), 164–172 (2012). [CrossRef]
5. C. Darne, Y. Lu, and E. M. Sevick-Muraca, “Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update,” Phys. Med. Biol. 59(1), R1–R64 (2014). [CrossRef]
6. A. D. Klose, B. J. Beattie, H. Dehghani, L. Vider, C. Le, V. Ponomarev, and R. Blasberg, “In vivo bioluminescence tomography with a blocking-off finite-difference SP3 method and MRI/CT coregistration,” Med. Phys. 37(1), 329–338 (2010). [CrossRef]
7. 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]
8. 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]
9. H. Yan, Y. Lin, W. C. Barber, M. B. Unlu, and G. Gulsen, “A gantry-based tri-modality system for bioluminescence tomography,” Rev. Sci. Instrum. 83(4), 043708 (2012). [CrossRef]
10. B. Li, F. Maafi, R. Berti, P. Pouliot, E. Rheaume, J. C. Tardif, and F. Lesage, “Hybrid FMT-MRI applied to in vivo atherosclerosis imaging,” Biomed. Opt. Express 5(5), 1664–1676 (2014). [CrossRef]
11. J. Zhang, D. F. Chen, J. M. Liang, H. D. Xue, J. Lei, Q. Wang, D. M. Chen, M. Meng, Z. Y. Jin, and J. Tian, “Incorporating MRI structural information into bioluminescence tomography: system, heterogeneous reconstruction and in vivo quantification,” Biomed. Opt. Express 5(6), 1861–1876 (2014). [CrossRef]
12. A. J. Chaudhari, A. A. Joshi, F. Darvas, and R. M. Leahy, “A method for atlas-based volumetric registration with surface constraints for optical bioluminescence tomography in small animal imaging,” in Proc. SPIE, (2007), p. 651024.
13. M. Baiker, J. Milles, J. Dijkstra, T. D. Henning, A. W. Weber, I. Que, E. L. Kaijzel, C. W. Lowik, J. H. Reiber, and B. P. Lelieveldt, “Atlas-based whole-body segmentation of mice from low-contrast Micro-CT data,” Med Image Anal 14(6), 723–737 (2010). [CrossRef]
14. M. Baiker, M. Staring, C. W. Lowik, J. H. Reiber, and B. P. Lelieveldt, “Automated registration of whole-body follow-up micro-CT data of mice,” Med. Image Comput. Assist. Interv. 14, 516–523 (2011). [CrossRef]
15. D. Xiao, D. Zahra, P. Bourgeat, P. Berghofer, O. A. Tamayo, H. Green, M. C. Gregoire, and O. Salvado, “Mouse whole-body organ mapping by non-rigid registration approach,” in Proc. SPIE, (2011), p. 79650E.
16. S. Ren, H. Hu, G. Li, X. Cao, S. Zhu, X. Chen, and J. Liang, “Multi-atlas registration and adaptive hexahedral voxel discretization for fast bioluminescence tomography,” Biomed. Opt. Express 7(4), 1549–1560 (2016). [CrossRef]
17. A. D. Klose and N. Paragas, “Automated quantification of bioluminescence images,” Nat. Commun. 9(1), 4262 (2018). [CrossRef]
18. H. Wang, D. B. Stout, and A. F. Chatziioannou, “Estimation of mouse organ locations through registration of a statistical mouse atlas with micro-CT images,” IEEE Trans. Med. Imaging 31(1), 88–102 (2012). [CrossRef]
19. H. Wang, D. B. Stout, and A. F. Chatziioannou, “A deformable atlas of the laboratory mouse,” Mol. Imaging Biol. 17(1), 18–28 (2015). [CrossRef]
20. B. Zhang, W. Yin, H. Liu, X. Cao, and H. Wang, “Bioluminescence tomography with structural information estimated via statistical mouse atlas registration,” Biomed. Opt. Express 9(8), 3544–3558 (2018). [CrossRef]
21. Kwangmoo Koh, Seung-Jean Kim, and S. Boyd, “l1_ls: Simple Matlab solver for l1-regularized least squares problems,” https://web.stanford.edu/~boyd/l1_ls/.
22. S. J. Kim, K. Koh, M. Lustig, S. Boyd, and D. Gorinevsky, “An interior-point method for large-scale l1-regularized least squares,” IEEE Journal on Selected Topics in Signal Processing 1(4), 606–617 (2007). [CrossRef]
23. M. E. Tipping, “Sparse bayesian learning and the relevance vector machine,” Journal of Machine Learning Research 1, 211–244 (2001).
24. D. Wipf and S. Nagarajan, “Iterative reweighted l1 and l2 methods for finding sparse solutions,” IEEE Journal of Selected Topics in Signal Processing 4(2), 317–329 (2010). [CrossRef]
25. Z. Zhang and B. D. Rao, “Sparse signal recovery with temporally correlated source vectors using sparse Bayesian learning,” IEEE Journal of Selected Topics in Signal Processing 5(5), 912–926 (2011). [CrossRef]
26. J. Feng, K. Jia, Z. Li, B. W. Pogue, M. Yang, and Y. Wang, “Bayesian sparse-based reconstruction in bioluminescence tomography improves localization accuracy and reduces computational time,” J. Biophotonics 11(4), e201700214 (2018). [CrossRef]
27. 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]
28. K. Liu, Y. J. Lu, J. Tian, C. H. Qin, X. Yang, S. P. Zhu, X. A. Yang, Q. S. Gao, and D. Han, “Evaluation of the simplified spherical harmonics approximation in bioluminescence tomography through heterogeneous mouse models,” Opt. Express 18(20), 20988–21002 (2010). [CrossRef]
29. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef]
30. M. Jermyn, H. Ghadyani, M. A. Mastanduno, W. Turner, S. C. Davis, H. Dehghani, and B. W. Pogue, “Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography,” J. Biomed. Opt. 18(8), 086007 (2013). [CrossRef]
31. 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]
32. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–61 (2013). [CrossRef]
33. Z. Zhang and B. D. Rao, “Extension of SBL algorithms for the recovery of block sparse signals with intra-block correlation,” IEEE Transactions on Signal Processing 61(8), 2009–2015 (2013). [CrossRef]
34. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007). [CrossRef]
35. J. Yu, B. Zhang, I. I. Iordachita, J. Reyes, Z. Lu, M. V. Brock, M. S. Patterson, J. W. Wong, and K. K. Wang, “Systematic study of target localization for bioluminescence tomography guided radiation therapy,” Med. Phys. 43(5), 2619–2629 (2016). [CrossRef]
36. J. Liu, Y. Wang, X. Qu, X. Li, X. Ma, R. Han, Z. Hu, X. Chen, D. Sun, R. Zhang, D. Chen, D. Chen, X. Chen, J. Liang, F. Cao, and J. Tian, “In vivo quantitative bioluminescence tomography using heterogeneous and homogeneous mouse models,” Opt. Express 18(12), 13102–13113 (2010). [CrossRef]
37. M. A. Naser, M. S. Patterson, and J. W. Wong, “Algorithm for localized adaptive diffuse optical tomography and its application in bioluminescence tomography,” Phys. Med. Biol. 59(8), 2089–2109 (2014). [CrossRef]
38. H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization,” Opt. Express 18(3), 1854–1871 (2010). [CrossRef]
39. X. He, J. Liang, X. Wang, J. Yu, X. Qu, X. Wang, Y. Hou, D. Chen, F. Liu, and J. Tian, “Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method,” Opt. Express 18(24), 24825–24841 (2010). [CrossRef]
40. Q. Zhang, X. Chen, X. Qu, J. Liang, and J. Tian, “Comparative studies of lp-regularization-based reconstruction algorithms for bioluminescence tomography,” Biomed. Opt. Express 3(11), 2916–2936 (2012). [CrossRef]
41. S. Jiang, J. Liu, Y. An, G. Zhang, J. Ye, Y. Mao, K. He, C. Chi, and J. Tian, “Novel L2,1-norm optimization method for fluorescence molecular tomography reconstruction,” Biomed. Opt. Express 7(6), 2342–2359 (2016). [CrossRef]
42. Y. Gao, K. Wang, Y. An, S. Jiang, H. Meng, and J. Tian, “Nonmodel-based bioluminescence tomography using a machine-learning reconstruction strategy,” Optica 5(11), 1451–1454 (2018). [CrossRef]
43. J. Yu, C. Dai, X. He, S. S. Guo H, and Y. Liu, “Bioluminescence tomography based on one-dimensional convolutional neural networks,” Frontiers in Oncology 11, 760689 (2021). [CrossRef]
