A Multilevel Probabilistic Cerenkov Luminescence Tomography Reconstruction Framework Based on Energy Distribution Density Region Scaling

Cerenkov luminescence tomography (CLT) is a promising non-invasive optical imaging method with three-dimensional semiquantitative in vivo imaging capability. However, CLT itself relies on Cerenkov radiation, a low-intensity radiation, making CLT reconstruction more challenging than other imaging modalities. In order to solve the ill-posed inverse problem of CLT imaging, some numerical optimization or regularization methods need to be applied. However, in commonly used methods for solving inverse problems, parameter selection significantly influences the results. Therefore, this paper proposed a probabilistic energy distribution density region scaling (P-EDDRS) framework. In this framework, multiple reconstruction iterations are performed, and the Cerenkov source distribution of each reconstruction is treated as random variables. According to the spatial energy distribution density, the new region of interest (ROI) is solved. The size of the region required for the next operation was determined dynamically by combining the intensity characteristics. In addition, each reconstruction source distribution is given a probability weight value, and the prior probability in the subsequent reconstruction is refreshed. Last, all the reconstruction source distributions are weighted with the corresponding probability weights to get the final Cerenkov source distribution. To evaluate the performance of the P-EDDRS framework in CLT, this article performed numerical simulation, in vivo pseudotumor model mouse experiment, and breast cancer mouse experiment. Experimental results show that this reconstruction framework has better positioning accuracy and shape recovery ability and can optimize the reconstruction effect of multiple algorithms on CLT.


INTRODUCTION
Cerenkov radiation (CR) is a physical phenomenon that occurs when a particle traveling through an object travels faster than the speed of light in the medium (1). The first radiographic imaging using CR was made in 2009 and is known as Cerenkov luminescence imaging (CLI) (2). Since its inception, CLI has been widely used in surgical guidance, drug development, endoscopic imaging, tumor detection, and other fields (3)(4)(5)(6)(7). The most significant advantage of CLI over other optical imaging methods is that it can use many approved radioactive sources for clinical imaging (8)(9)(10). However, CLI is a planar imaging method and cannot obtain the depth information and three-dimensional (3D) distribution of radioactive sources. Therefore, a new optical imaging method, Cerenkov luminescence tomography (CLT), combined with CLI and 3D anatomical imaging modality, has been developed. Compared with CLI, CLT can obtain the internal and external contour or boundary of biological tissues with the 3D anatomical imaging modality and determine the 3D spatial distribution of radioactive source in biological tissues (11)(12)(13)(14)(15).
However, the Cerenkov photon energy mainly concentrated in the short-wavelength band with a high scattering in biological tissues, which leads to the difficulty for CLT reconstruction (16). Therefore, compared with bioluminescence tomography (BLT) and fluorescence molecular tomography (FMT), CLT requires more prior information constraints and more robust algorithms to optimize its solution (17,18). Therefore, researchers have carried out a series of work in algorithms and feasible region constraints for CLT reconstruction.
From the algorithm's point of view, improving the current reconstruction algorithm is a research focus, such as L 1 norm regularization (19), L 2 norm regularization (11), LP (0<p<1) norm regularization, and other regularization methods (12), that has been used in the CLT field. Although these regularization methods can be applied to CLT, they depend on the selection of regularization parameters. In addition to regularization algorithms, there are other areas of inverse problem algorithms; non-regularization methods such as orthogonal matching pursuit (OMP) (20) do not require regularization parameter selection, but they have poor performance in CLT reconstruction.
From the perspective of the feasible region, some feasible region iterative methods, such as iterative shrinking permissible region (ISPR) (21), three-way decision (TWD) (22), and feature extraction from the autoencoder (3), have already been applied to optical molecular imaging. These methods constrain the solution space by regions of interest (ROIs) to obtain better reconstruction results. However, the current feasible region method has the following shortcomings: first, the feasible region only shrinks and may lose some nodes; second, the final reconstruction results are only affected by the previous one, ignoring the intermediate process; and third, errors in the iterative process will be passed to the subsequent iterative methods, leading to polluted results. In general, most of the reconstruction algorithms or strategies are limited by the weak signal intensity of CR, which leads to the difficulty of solving the inverse problem, or the parameter tuning of the reconstruction algorithm itself needs to be carried out at a high cost. From what has been discussed above, many inverse problem algorithms or feasible region strategies are challenging to apply to the CLT field.
To optimize the feasible region method in CLT and improve the performance of traditional algorithms in CLT reconstruction, this paper proposed a multilevel probabilistic energy distribution density region scaling (P-EDDRS) framework for CLT. In this framework, L 2 norm error rate and cosine similarity were used to evaluate the quality of each iteration reconstruction result, and normalized weight was assigned to each reconstruction result according to the evaluation. This normalized weight represents the probability that the corresponding iteration result is the final result. By this weight, the initial iteration reconstruction result becomes a probabilistic result. Then, the probabilistic result is regarded as random variables distributed in 3D space, and ROI for the next iteration is divided according to the distribution density of these random variables. Besides, to stabilize the rate of ROI change, the formula proposed by Naser et al. was introduced (23). After several iterations, all the initial reconstruction source distribution and the corresponding normalized weight are weighted to get the final Cerenkov source distribution.
To evaluate the performance of the P-EDDRS frame in CLT, several groups of numerical simulations and in vivo experiments were implemented. The optimized Lasso and Least Square QRfactorization (LassoLSQR) algorithm based on L 1 norm regularization (24), Tikhonov regularization algorithm based on L 2 norm regularization (25), damped singular value decomposition (DSVD) algorithm based on SVD (26), and OMP algorithm based on matching pursuit (20)-these algorithms include regularization algorithm and greedy algorithm-was used to evaluate the performance of the framework combined with various reconstruction algorithms. In addition, to compare this CLT framework's performance with those of other ROI methods, ISPR and TWD methods are introduced (21,22). The parameters of all reconstruction algorithms and methods take their default values. The results prove that the P-EDDRS framework for CLT can combine various algorithms to reconstruct radioactive sources of different sizes and shapes, which has higher applicability. In addition, in vivo experiments show that the framework is still reliable and stable in living animals.

Inverse Problem
Most of the energy of CR is concentrated in the short wavelength band, characterized by high scattering and low absorption in biological tissues. So the diffusion approximation (DA) model can be used as the mathematical basis to describe the Cerenkov photon transport process in tissues. DA models with Robin boundary conditions are often described as follows (27): where F(r) denotes the Cerenkov photon flow rate at the point r in the region W, m a and m s are the absorption coefficient and scattering coefficient of a tissue, g is anisotropy factor, and D(r) denotes the diffusion coefficient at position r. The symbol ∇ is used for the differential operator of vector, ?W is the set of the surface (boundary) points, and x is the point on the surface of a tissue. R f is the internal indicator of refraction of the tissue, and n is the unit normal vector whose direction is from the inside of the biological tissues to the outside of ?W. Furthermore, the continuous space in biological tissues can be discretized into finite units by using the finite element method (FEM). By using FEM in solving Eq. (1), a reduced linear relationship between the unknown Cerenkov source distribution in the tissue and the surface photon flow rate can be obtained: where A is the CLT system matrix, and it gives the mathematical process of Cerenkov photon propagation in the tissue, B represents the Cerenkov photon flow rate vector on the surface of biological tissues measured by a susceptible CCD camera, length M of B represents the semaphore measured, X represents the unknown Cerenkov source distribution vector in the organism, and length N of X represents number of grid nodes of FEM.
In essence, most of feasible region iterative reconstruction methods are essentially subtracted columns from A to simplify the process of solving Eq. (2). For example, the TWD method divides the initial reconstruction region into positive domain (POS), negative domain (NEG), and boundary domain (BND); combines POS and BND into ROI; and deletes the system matrix column corresponding to NEG. In the ISPR method, each iteration's reconstruction results are arranged in descending order according to the energy intensity of nodes, and the columns of system matrix corresponding to nodes with lower energy are removed. Those iterative reconstruction methods based on feasible region have defects as mentioned in the Introduction. To optimize the feasible region method in CLT and improve the performance of traditional algorithms in CLT reconstruction, the P-EDDRS framework is proposed in this paper.

Probabilistic Energy Distribution Density Region Scaling Framework for Cerenkov Luminescence Tomography
The P-EDDRS framework for CLT mainly includes the following steps: 1. Reconstruct the CR source based on ROI (initial ROI is global).
2. Evaluate the quality of each iteration reconstruction result and generate probabilistic source distributions based on the evaluation.
3. According to the reconstruction source distribution density, determine the subsequent ROI.
4. Dynamically change the next ROI size based on the energy intensity of the current result.
5. Evaluate whether the cutoff conditions are met and continue if so; otherwise, return to step 1.
6. The results are post-processed and weighted to get the final Cerenkov source distribution when the iteration is complete. Figure 1 shows the flowchart of the P-EDDRS framework for CLT.

Evaluate Reconstruction Result and Generate Probabilistic Source Distribution
To evaluate the error of each CLT reconstruction radiation source distribution, this section performed the following actions.
First, the index named Index_Sp is initialized based on the global number of nodes. Next, according to Index_Sp, the reconstruction algorithm mentioned above is used to reconstruct the current region (the first time is the whole FIGURE 1 | Flowchart of P-EDDRS framework for CLT. P-EDDRS, probabilistic energy distribution density region scaling; CLT, Cerenkov luminescence tomography. region) radiation source distribution, and the first iteration reconstruction source distribution X 1 is obtained. According to Eqs (3) and (4), the L 2 norm error rate E L2(x 1 ) and the cosine similarity E Cos(x 1 ) of X 1 can be obtained, respectively.
where i is the number of iterations between 1 and the maximum number of iterations L max 50. For now, evaluation of the reconstruction is obtained through Eqs (3) and (4). In order to assign the corresponding weight to each reconstruction source distribution in the overall framework, the weight value of each iteration reconstruction source distribution is introduced. The weight here essentially represents the probability value that the result of each iteration is the final result, so we can call this weight value as the probability weight value. Furthermore, from the reconstruction perspective, the L 2 norm error is an inversely proportional evaluation. To achieve this effect, introduce Eq. (5).
The other indicator E cos is different from E L2 . The closer the similarity is to 1, the more reliable the reconstruction source distribution is. It is a proportional relationship, which can be represented by Eq. (6).
Now, two metrics can use to evaluate each outcome. To integrate these two indexes, Eq. (7) can be introduced. Eq. (7) gives the result of each reconstruction source distribution corresponding probability weight value. When all iterations are completed, each iteration's initial reconstruction source distribution is multiplied by the corresponding probability weight, and the sum is the final reconstruction source distribution.
Therefore, the normalized probabilistic result of each reconstruction source distribution can be expressed by Eq. (8). Through Eq. (8), the initial reconstruction source distribution is transformed into probabilistic reconstruction source distribution with weights. The distribution can be combined with the corresponding grid coordinates as a group of random variables distributed in 3D space to determine the ROI region in the next section.

Determine the Region of Interest
After the probabilistic source distribution of one reconstruction in Section 2.3 is obtained, the feasible area of the subsequent reconstruction needs to be determined. We assume that the probabilistic source distribution X iP is a set of random variables whose mathematical expectation is m and variance is s 2 . For any positive number K, Chebyshev's inequality holds as in Eq. (9).
Further, K is expressed as different values: It can be seen from the above formula that the distribution of random variables has its inherent trend. According to this feature, assume that all nodes with corresponding grid coordinates are random variables distributed in 3D space. The distribution of random variables in 3D space also has the rule of Eq. (10). In other words, the distribution of all reconstruction source distribution tends to be close to its mathematical expectation, which is the basis for determining ROI. According to the characteristics of Eq. (10), cuboid can be used as the form of ROI. Now, we have a center point, the length of the sides in the three directions and the deflection angle concerning the coordinate axes are determined, and the form of ROI (cuboid) in space can be determined. The average probability of nodes can be obtained by Eq. (11) as the center of the ROI region.
where N represents the total number of nodes in the current iteration; x t , y t , and z t represent the coordinates of the point t; x iP(t) is the tth element in X iP ; x P , y P , and z P represent the probabilistic average value of the corresponding coordinate. Once the center is determined, the deflection angle and side length of ROI need to be determined. The covariance matrix M cov is introduced to solve the side length and deflection angle of ROI.
where s 2 x , s 2 y , and s 2 z represent the probability variance of each coordinate direction; and Cov xy , Cov yz , and Cov xz represent the covariance between different coordinate directions. These two variables are represented by Eqs (13) and Eq. (14): All parameters are multiplied by N/(N − 1) to ensure unbiased estimates of the variance and covariance. Through the above calculation, the construction of the overall probability result covariance matrix is completed. Further, the eigenvalues Val and eigenvectors Vec of the matrix can be obtained. Val and Vec are in the form of matrices, in which the three diagonal elements of Val represent the corresponding eigenvalues in the three coordinate directions. At the same time, the Vec is the eigenvector in the three coordinate directions. To realize dynamic scaling of the region, the ROI side length R can be defined as in Eq. (15): where Size is the coefficient used to control the ROI, which is initially set to 1 and will be refreshed in Section 2.5. By Eqs (11) negative signs are determined, the vector direction from the origin to each vertex can be deflected by multiplying the corresponding feature vectors. The center point obtained by Eq. (11) can be added to obtain the eight vertices of the current ROI. At the same time, we can dynamically adjust the scaling rate of next ROI by multiplying the size variable introduced in Section 2.5 with Val.
When the ROI is determined, nodes among this ROI can be detected, and calculate the spatial distance between the nodes in the ROI region and the ROI center and arrange them in ascending order to get the ROI node index named Index_Sp_ROI. Figure 2 shows how the ROI is generated. Figures 2A-C show examples of one of the ROI division processes, and Figure 2D shows examples of all ROI tracks of the overall process.

Dynamically Change the Region of Interest Size
The previous step identified the ROI area, but the ROI rate of change will be significantly accelerated during the mid to late iteration phases. This phenomenon will cause some nodes to be directly classified outside the ROI, which will lead to the loss of the morphological information of the radiation source. To better retain the nodes that may be the radiation source and more information about the energy of radiation source, the change rate of ROI should be controlled. So this paper introduces the formula of Naser et al., as follows (23) where Num f is the expected number of nodes remaining in the final iteration, because the grid is used as the division unit of spatial structure, and this value is set as 4 (a tetrahedron).
Cut_Num is size of Index_Sp_ROI at first iteration, and ß represents the attenuation coefficient. Now, a judgment condition between the size relationship between Cut_Num/ß 2 and the Cut_Num has been set: if the quotient is greater than 2, the Size in Eq. (15) is set to 2; if it is less than 1, the Size is set to 0.5; in other cases, the Size is set to 1. This variable enables the ROI to shrink or expand.
After the change of ROI region is completed, the current Index_Sp is arranged in descending order according to the energy intensity value of the CLT reconstruction source distribution. The new Cut_Num is made by dividing the current Cut_Num by ß and rounding the number toward positive.
The previous Cut_Num nodes of Index_Sp are taken as the index of the current threshold shrinking to generate Index_Sp_Descend.
Concatenate index Index_Sp_Descend after Index_Sp_ROI, and a new index is obtained, which contains the nodes of the ROI region and the nodes with the threshold shrunk. After the duplicate elements of this new index are removed, use Cut_Num to cut off the new index, and the remaining elements are the Index_Sp, which will be used for the following CLT reconstruction. When the remaining elements are still greater than 1 and the number of iterations does not reach L max , return to Section 2.3.

Post-Processing and Acquisition of the Final Result
At the end of the overall iteration, all iteration reconstruction results to X i , and the corresponding probability error rate P Err(Xi) has been obtained. Although all reconstruction source distribution has been recorded, some iterations have significantly higher probabilistic error rates. From the description by Ding et al., the iteration result error evaluation of 3D reconstruction tends to be in the form of Gaussian distribution (13). Here, the error evaluation can be used as the criterion for evaluating the results, so the mean value and standard deviation for E L2(Xi) and E cos(Xi) can be calculated. According to the definition of Gaussian distribution, only the iteration source distribution that is centered on the mean and within one times the standard deviation is retained. The intersection of the two Gaussian filtering results is the result of some iterations that need to retain. Once this step is complete, use Eqs (5)-(7) again to generate the new parameters and P' Err(Xi) . The final CLT reconstruction source distribution X final can be obtained by Eq. (17).
where K is the number of remaining results after the completion of filtering treatment in this section.

EXPERIMENTS AND RESULTS
The CLT reconstruction experiments were performed in MATLAB 2020B and run on a desktop computer with a 3.00-GHz Intel Core i5 CPU and 32 GB of memory.
To verify and systematically evaluate the performance and characteristics of the CLT reconstruction framework proposed in this paper, several groups of CLT numerical simulation experiments and CLT in vivo experiments were designed. The location error (E L ), the dice coefficient (Dice), the tetrahedral volume ratio (R V ), and the global relative residual (R R ) are introduced as quantitative evaluation indicators.
The location error E L is defined as the Euclidean distance between the reconstruction source distribution center point coordinates (x,y,z) and the actual radiation source coordinates (x 0' y 0' z 0 ): Dice coefficient is used to evaluate the degree of shape similarity (overlap) between the reconstructed source distribution area R and the actual radiation source distribution area T. The closer it is to 1, the higher the similarity between R and T is. The following formula defines Dice coefficient: The tetrahedral volume ratio R V is defined as the ratio of the tetrahedral volume V T of the actual radiation source distribution to the tetrahedral volume V R of the reconstruction source distribution. Same as Dice, when R V is closer to 1, the higher the size similarity of the radiation source. The following formula calculates R V : The final theoretical error of the evaluation result of global relative residual R R is calculated by the following equation: the lesser the R R , the smaller the theoretical error.

Numerical Simulations
In this section, a non-homogeneous digital mouse model show in Figure 3 is applied to verify the performance of the CLT reconstruction framework (28). To reduce the computational complexity and performance cost, the head and tail of the digital rat were removed, and only the trunk containing the main organs was retained, which was composed of the heart, lung, liver, stomach, and kidney, while the rest were muscle tissues. 2-[ 18 F]-Fluoro-2-deoxy-D-glucose ( 18 F-FDG) was the Cerenkov radioactive source used. In the CLT numerical simulation, Cerenkov photons generated by 18 F radioactive source in tissues were simulated by GEANT4, and the transmission process of these Cerenkov photons in tissues was simulated by MOSE (29)(30)(31)(32).
The CR signal with wavelength of 630 nm was collected for reconstruction. The optical parameters of each organ and tissue are given in reference (33), as shown in Table 1. In the CLT numerical simulation experiment, the radiation source is distributed in two positions, namely, (9.5, 15.5, 25) mm and (20,8,15) mm, using different shapes and sizes. One set of shapes is shown in Figures 3A-D show the surface Cerenkov photon energy distribution maps corresponding to Figures 3A, B.
Four groups of experiments were designed in this section. The first group of experiments uses Tikhonov, DSVD, LassoLSQR, and OMP algorithms mentioned above to reverse reconstruct the radiation source in combination with the framework to verify the feasibility of this framework. The second set of experiments compared this framework with ISPR and TWD iterative reconstruction methods to verify the framework's effectiveness in CLT. In the third and fourth groups, radiation sources of different sizes and shapes were selected for reconstruction to verify the stability and robustness of the framework in CLT.

The Experiment of Feasibility Verification
This experiment was designed to verify the feasibility of the proposed framework in CLT and the versatility of different algorithms. The digital mouse has been discretized into a grid of 12,831 nodes and 66,085 tetrahedrons using Amira (Visage Imaging, Australia). A spherical radiation source with a radius of 0.8 mm was placed at coordinates (9.5, 15.5, 25) mm as shown in Figure 3A. The results are shown in Figure 4, and the quantitative indexes are shown in Table 2.
In Figure 4, the radiation source is marked in red, and the reconstruction source distribution is shown in a purple grid in all 3D views. In all axial views, the actual location of the radiation source is marked by a red arrow. It can be seen that in all four algorithms, the reconstruction source distribution is significantly improved after the framework is used, and different algorithms can realize convergence after several iterations and obtain better reconstruction results. Specifically, Tikhonov and DSVD algorithms have poor convergence in one-step reconstruction results and cannot determine the actual radiation source position. With P-EDDRS, the correct radiation source location can be obtained. LassoLSQR and OMP algorithms' initial result is incorrect and difficult to converge into a single region. After using the framework, the reconstruction results can be corrected. In addition, in terms of the degree of shape recovery, this framework can better restore the shape of the radiation source, which can be proved from the pictures and the Dice coefficient. It can also be seen from RR in Table 2 that the results of this framework are closer to the theoretically more correct results.

The Experiment of Efficiency Verification
This experiment was designed to compare the performance of the proposed framework with two other typical ROI shrinking methods: TWD and ISPR (21,22) in CLT reconstruction. The grid used and the radiation source placement are the same as those of the last experiment. The Tikhonov and DSVD algorithms mentioned above were used to reverse reconstruct the radiation source combined with the framework. The results are shown in Figure 5, and the quantitative indexes are shown in Table 3.
In Figure 5, the radiation source is marked in red, and the reconstruction source distribution is shown in a purple grid in all 3D views. In all axial views, the actual location of the radiation source is marked by a red arrow. As shown in Figure 5A, ISPR leads Tikhonov algorithm results to tend to the point with high energy intensity value, and the correct radiation source depth cannot be obtained. At the same time, in Figure 5B, similar trends tend to be of high energy intensity but ignore the shape. Because ISPR determines the ROI region based on node energy drop ranking, spatial information between nodes is ignored, resulting in the follow-up's discontinuous morphological  distribution, leading to multiple high energy points. The value of R R under TWD+Tikhonov condition is smaller than that under P-EDDRS+DSVD condition. However, the latter result is more accurate whether evaluated from other quantitative indicators, 3D views, or cross-section views. As shown in Figures 5C, D, TWD can obtain the actual radiation source depth relatively closely. However, the results of TWD depend entirely on the previous iteration, which will lead to the accumulation of error results and finally get the error distribution. In Figures 5E, F, our framework can achieve better results among the three methods in terms of the results of the two algorithms, whether it is positioning error, Dice coefficient representing morphology, global relative residual, or the volume of reconstructed radiation source distribution. In addition, the depth and distribution of the radiation source are closer to the actual location.

The Stability and Robustness Experiment I
This experiment is to verify the stability and robustness of the P-EDDRS framework for different radiation source sizes in CLT. Four spherical radiation sources with different radii and sizes are used as targets. The digital mouse has been discretized into a grid of 28,463 nodes and 159,957 tetrahedrons in the same way. The spherical radiation source was placed at coordinates (20,8,15) mm as shown in Figure 3B. Different from the first two, four radiation sources of different sizes have been placed at grid in this  Figure 6, and the quantitative indexes are shown in Table 4.
In Figure 6, the radiation source is marked in red, and the reconstruction source distribution is shown in a purple grid in all 3D views. In all axial views, the actual location of the radiation source is marked by a red arrow. As shown from the results in Figure 6, our framework can achieve better performance in four different sizes of experiments, and all the indicators are satisfactory. From all the results, the E L of all results is between 0.50 and 0.77, and the Dice coefficient remained above 71%. All but the minor radiation source is close to the actual volume. Because the finite element grid method is adopted in this paper, it has a specific size. If the radiation source is too small, the grid cannot fit the radiation  source volume well. This phenomenon leads to some errors in the reconstruction of small size radiation sources. This is the limitation of the grid reconstruction method. This experiment proves that different reconstruction algorithms can estimate the size of radiation sources effectively under this framework in CLT.

The Stability and Robustness Experiment II
This experiment is similar to the last in verifying the stability and robustness of the P-EDDRS framework for different radiation source shapes in CLT. Four different shapes of radiation sources were placed in the grid: a sphere with a radius of 1.25, a cube with sides of 2.5, a cylinder with a radius of 1.25 and a height of 2.5, and an ellipsoid with a = 2 and b = c = 1.25. The grid is the same as in the last experiment. The radiation source was placed at coordinates (20,8,15) mm. The results are shown in Figure 7, and the quantitative indexes are shown in Table 5.
In Figure 7, the radiation source is marked in red, and the reconstruction source distribution is shown in a purple grid in all 3D views. In all axial views, the actual location of the radiation source is marked by a red arrow. From all the results, the E L of all results is between 0.50 and 0.86, and the Dice coefficient remained above 65%. All of the reconstruction source distribution volume is close to the actual radiation source volume. This experiment proves that different reconstruction algorithms can estimate the shape of radiation sources effectively under this framework in CLT.

In Vivo Experiments
To verify the feasibility and performance of our framework in an actual CLT situation, two groups of in vivo experiments are designed. The first group verifies the reconstruction ability of the framework in the case of relatively deep radioactive sources by using the pseudotumor model, and the second group verifies the reconstruction ability of the framework in the case of shallow radioactive sources by using subcutaneous breast cancer. Two approximately 7-week-old female nude mice (BALB/c Nude) served as the imaging model. All animal procedures were performed under isoflurane gas anesthesia (2% isoflurane-air mixture) to minimize the suffering of the mice. Animal experiments comply with the Regulations on the Management of Experimental Animals. All procedures follow the Animal Ethics Committee of the Northwest University of China (No. NWU-AWC-20210901M). The optical data were acquired using the iXon Ultra electron double CCD camera manufactured by ANDOR (Northern Ireland). The X-ray source is the L9181-02 microfocus ray source, and the X-ray detector is C7942CA-22, all manufactured by HAMAMATSU (Japan). The optical lens is EF 24 mm f/1.4L II USM manufactured by Canon (Japan). The band-pass filter is FF01-630/92-25 manufactured by Semrock (USA). The CCD camera is pre-cooled to −80°C, and the fuzzy local information C-means and curvature-driven diffusion (FLICMCDD) method in reference (34) is used to reduce the influence of noise signals. Exposure time is set to 5 min, the gain value is set to 300, shift speed is set to 13 ms, and the read rate is set to 1 MHz at 16 bits. In in vivo experiments, 18 F-FDG is used as the Cerenkov radioactive source. The data acquisition equipment is shown in Figure 8A.

The Experiment of Pseudotumor
The pseudotumor model used in this study was made of plastic and was a cylindrical tube with the most extended length of 2.3 mm, the widest width of 1.1 mm, and the highest length of 1.5 mm, as shown in Figure 8B. This pseudotumor model was injected with about 800 ± 50 mCi of 18 F-FDG as Cerenkov radioactive source. The pseudotumor model was implanted in the mouse abdominal cavity, at the back of the left inner lobe of the liver, close to the liver, and implantation depth is deeper. The mouse was dispersed into a grid of 13,475 nodes and 69,923 tetrahedrons using FEM as shown in Figure 8B. This grid also removed the head and tail to reduce computational complexity, leaving only the main organs in the trunk. The difference from the numerical simulation is that this grid contains bones. The corresponding optical parameters of animal tissues were the same with numerical simulations. In this experiment, the approximate central coordinates of the pseudotumor model are (19.68, 18.40, 19.20) mm. The Tikhonov and DSVD algorithms mentioned above were used to reverse reconstruct it in combination with the P-EDDRS framework. The reconstruction results are shown in Figure 8, and the quantitative indexes are shown in Table 6.
In Figure 8, the radiation source is marked in red, and the reconstruction source distribution is shown in a purple grid in all 3D views. In all axial views, the actual location of the radiation source is marked by a red arrow. As shown in Figures 8C, D, Tikhonov and DSVD algorithms in one-step reconstruction cannot obtain the correct radiation source depth in the in vivo CLT experiment. In addition, it can be seen from Figures 8C, D that Tikhonov and DSVD algorithms identify multiple radiation source distribution. However, in Figures 8E, F, the P-EDDRS framework can ensure that the distribution of the reconstructed radioactive sources tends to be a whole, and the center position of the reconstructed radioactive sources with the P-EDDRS framework was relatively close to the actual pseudotumor model. According to the quantitative results in      Table 6, the E L of all results is between 0.78 and 0.89, and the Dice coefficient directly rises from 0 in the initial result to more than 50%.

The Experiment of Subcutaneous Breast Cancer
In this experiment, a mouse implanted with 4T1 breast cancer cells was used as the target animal. The optical parameters, equipment, and experiment setting are the same as in Section 3.2.1. Special thanks to the Institute of Automation, Chinese Academy of Sciences for the data provided. The mouse was dispersed into a grid of 14,289 nodes and 71,838 tetrahedrons using FEM. The difference is that this experiment used the subcutaneous tumor and is not influenced by other tissues and organs. Therefore, to reduce the computational complexity, the tissues and organs inside the mouse were removed, and the mouse was regarded as a homogeneous structure composed of muscle tissues; 1 × 10 6 4T1 cells were subcutaneously injected into the back of the mouse. After 6 days of culture, about 800 ± 50 mCi of 18 F-FDG was injected through the tail vein as a radioactive source. After 40 min, a CCD camera was used, the exposure time was 5 min, and the image in Figure 9A was collected. The Tikhonov and DSVD algorithms mentioned above were used to reverse reconstruct it in combination with the framework. In this experiment, the approximate central coordinates of the tumor are (22.7, 24.4, 11.0) mm. Breast cancer tumor was used as the focus in this experiment, and PET is currently the gold standard for imaging tumors, so PET data were introduced as a reference in Figure 9A.
The figure of collected data, approximate tumor location obtained by PET, and reconstruction results are shown in Figure 9, and quantitative indicators are shown in Table 7.
In Figure 9, the radiation source is marked in red, and the reconstruction source distribution is shown in a purple grid in all 3D views. In all axial views, the actual location of the radiation source is marked by a red arrow. PET scan results were compared in this experiment to better compare with the actual situation, as shown in Figure 9A. It can be seen from Figure 9A that the noise signal collected in CLT contains many high-energy particles and thermal noise caused by prolonged exposure. This makes noise reduction of CLT particularly important and necessary. Therefore, this paper introduces the FLICMCDD denoising algorithm to preprocess the collected signals (34).
It can be seen from Figure 9B that the radiation source distribution obtained by Tikhonov algorithm in a single calculation is too sparse, with only a few nodes. The same conclusion can be reached in the volume ratio of quantitative results. Compared with Tikhonov algorithm, the radiation source distribution obtained by DSVD algorithm is larger than the actual tumor, because the noise obtained by this group of data is too large, resulting in severe artifacts. The comparison of Figures 9B-E shows better results. According to the quantitative results, Tikhonov and DSVD algorithms have improved radiation source positioning accuracy, size, and shape using this paper's framework in CLT. From all the results, the E L of all the results stays below 1. However, artifacts appeared in the results of different methods. This is because the grid of our group is too simplified and ignores the influence of other tissues of mouse on CR, which ultimately leads to a more significant distortion of the collected CR signal. This phenomenon also leads to a low Dice coefficient compared with numerical simulation experiments. However, it can be seen from the comparison of Figures 9D, E that the framework has a specific elimination effect on artifacts. This phenomenon happens because the reconstruction results with artifacts have a sizeable error rate, which leads to a significant decrease in their weight. In the final result, the intensity will be weakened accordingly, reducing the impact on the overall result, which will be shown as weak artifacts from the cross-section diagram. On the whole, the reconstructed radioactive source distribution in Figures 9C, E can be consistent with PET.

DISCUSSION AND CONCLUSIONS
This paper proposed a multilevel ROI-scaled CLT reconstruction framework based on probabilistic energy density. In this framework, all nodes with energy are regarded as random variables of 3D spatial distribution, and the cuboid ROI region is used to ensure spatial continuity between nodes. Optimize the defect of neglecting the spatial distribution of nodes in the iterative reconstruction method only based on energy intensity. In order to ensure a stable ROI regional change rate, a shrinking formula based on energy intensity is introduced to ensure that energy information and spatial information of nodes are considered simultaneously. By the idea of iterative reconstruction, the corresponding probability value is assigned to the result of each reconstruction and dynamically refreshed in each iteration to avoid the contamination of the global result by some iteration errors such as radioactive source is too sparse and scattered or the artifact is too prominent in CLT.
To verify and evaluate the feasibility and performance of the proposed framework in CLT, numerical simulation and in vivo experiments are introduced. Four reconstruction algorithms, Tikhonov, DSVD, LassoLSQR, and OMP, were used for qualitative and quantitative analyses and comparison. The following conclusions can be obtained from the experimental results: first, the feasibility experiment shows that our framework can improve the radioactive source positioning accuracy of different algorithms in CLT; all E L values stay below 1; second, our frame has advantages in positioning accuracy and radioactive source shape recovery ability for CLT compared with other one-step reconstruction or iterative reconstruction method, the Dice coefficient of all numerical simulation experiments is above 0.65 and above 0.43 in vivo experiments; third, the experimental results of different sizes and shapes of radioactive sources show that our frame is robust and can better obtain the size and shape of radioactive sources for CLT; finally, in vivo experiments for CLT verified the feasibility of this framework in the detection of the radioactive source in live animals. It is worth noting that the reconstruction result of in vivo experiment is slightly worse than that of the numerical simulation experiment. This is because in in vivo experiments, the CCD camera collected a relatively weak CR signal and a large number of high-energy rays produced by radionuclides and thermal noise of the CCD camera itself. Furthermore, the deviation will inevitably occur in the process of energy mapping to the grid.
There are still some deficiencies in the framework. First of all, this framework is based on the idea of iterative multiple times. Compared with a one-step reconstruction, the calculation cost is high. Second, this framework is based on the finite element mesh method, resulting in inconsistent reconstruction performance for different shapes and sizes of radioactive sources. In addition, the cube ROI does not fit the irregular radioactive source well, leading to the low Dice coefficient in some cases. Subsequent work will attempt to reconstruct using other structures such as voxels or point clouds and try ROI of other shapes to overcome current deficiencies.    On the whole, the P-EDDRS frame treats the nodes resulting from each reconstruction as 3D random variables. According to the distribution characteristics of random variables, cuboid ROI was used to delimit molecular regions to ensure the spatial continuity of reconstruction results. Shrinking formula based on energy intensity is introduced to ensure that energy node information and spatial information jointly control regional ROI changes. Dynamic probabilistic results can guarantee the correctness of reconstruction results. The reconstruction results of numerical simulations and in vivo experiments demonstrate that our framework can improve the ability of multiple reconstruction algorithms to locate the radioactive source and recover the shape of the radioactive source.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Ethics Committee of the Northwest University of China (No. NWU-AWC-20210901M).

AUTHOR CONTRIBUTIONS
XW: conceptualization, methodology, software, programing, validation, visualization, data curation, writing-original draft, and writing-review and editing. HG: conceptualization, writingreview and editing, funding acquisition, and data curation. JY: writing-review and editing, and funding acquisition. XLH: conceptualization and writing-review and editing. HY: funding acquisition. YH: project administration and funding acquisition. XWH: writing-review and editing, supervision, project administration, and funding acquisition. All authors contributed to the article and approved the submitted version.