GCR-Net: 3D Graph convolution-based residual network for robust reconstruction in cerenkov luminescence tomography

Cerenkov Luminescence Tomography (CLT) is a novel and potential imaging modality which can display the three-dimensional distribution of radioactive probes. However, due to severe ill-posed inverse problem, obtaining accurate reconstruction results is still a challenge for traditional model-based methods. The recently emerged deep learning-based methods can directly learn the mapping relation between the surface photon intensity and the distribution of the radioactive source, which e®ectively improves the performance of CLT reconstruction. However, the previously proposed deep learning-based methods cannot work well when the order of input is dis-arranged. In this paper, a novel 3D graph convolution-based residual network, GCR-Net, is proposed, which can obtain a robust and accurate reconstruction result from the photon intensity of the surface. Additionally, it is proved that the network is insensitive to the order of input. The performance of this method was evaluated with numerical simulations and in vivo experiments. The results demonstrated that compared with the existing methods, the proposed method can achieve e±cient and accurate reconstruction in localization and shape recovery by utilizing three-dimensional information.


Introduction
Cerenkov luminescence imaging (CLI) is an optical molecular imaging technology based on Cerenkov radiation, which is electromagnetic radiation emitted when a charged particle travels faster than light in the medium. 1,2Like other optical molecular imaging technologies, CLI has the essential advantage of high sensitivity and low cost.][5][6][7][8][9][10][11][12][13][14][15] At present, CLI has been used for tumor imaging, therapy monitoring, etc.However, CLI only performs planar imaging and cannot obtain the three-dimensional distribution of radionuclides, which hinders its further clinical application.
To trace the location and recover the distribution of radionuclides, the tomography method for CLI was proposed and named Cerenkov Luminescence Tomography (CLT). 16However, the reconstruction of CLT is a highly ill-posed reconstruction problem that makes it di±cult to obtain accurate reconstruction results.1][22][23][24] For example, Feng et al. took permissible region as prior information and proposed an optimal permissible region selection strategy. 257][28][29][30] Zhong et al. proposed the whole body reconstruction algorithm along a regularization path. 31Besides the strategies mentioned above, taking multispectral information into account was also proposed. 32,33Spinelli et al. proposed a multispectral approach that used distinctive information content at di®erent wavelengths in the reconstruction process, and the result demonstrated that multispectral information can improve the e±ciency of CLT. 34he recent success in deep learning demonstrated its great potential in optical imaging tomography.By employing a deep neural network, the mapping relation between the surface photon intensity and internal radioactive source can be learned directly.Zhang et al. ¯rst applied a multilayer fully connected neural network (MFCNN) to reconstruct the distribution of the radioactive source. 35Cao et al. proposed a CLT reconstruction framework based on a stacked denoising auto-encoder, in which the residual error was learned between the preliminary reconstruction result and the accurate result. 36hang et al. applied an attention mechanism-based locally connected network for CLT, which improved the accuracy and stability of the CLT reconstruction. 37For the reconstruction in bioluminescence tomography and °uorescence molecular tomography, Meng et al. proposed a locally connected network based on the K-Nearest Neighbor method, which achieved accurate morphological FMT reconstruction in a short time. 38Yu et al. proposed a one-dimensional convolutional neural network, which signi¯cantly reduced the number of parameters. 39The methods mentioned above take the node intensity as the input of the network, and the nodes are generated by ¯nite element method (FEM) software. 20,26However, the trained network is sensitive to the order of input, in other words, even for the same imaging model, if the order of input is disarranged, the trained network cannot work well.
In this paper, a 3D graph convolution-based neural network with a residual architecture for CLT reconstruction, named GCR-Net, is proposed to overcome the shortcomings of existing reconstruction methods based on deep learning.GCR-Net is divided into several residual blocks and each block consists of basic residual units.First, a threedimensional graph convolution operator is introduced to extract and aggregate features.Second, coordinates and the surface photon intensity are performed as the input, which suggests that threedimensional information is utilized in our network.It is noteworthy that the proposed network does not rely on a speci¯c order for input points because of the characteristic of the three-dimensional graph convolution operator, which means that an accurate result can be obtained even if the order of input is disarranged.
The rest of this paper is structured as follows: In Sec. 2, basic knowledge of CLT reconstruction and the details of GCR-Net are introduced.In Sec. 3, the numerical simulation experiments as well as the in vivo studies are illustrated to prove the e®ectiveness of our network.In Sec. 4, the discussion and conclusion are illustrated.

Methodology 2.1. Model-based CLT reconstruction
For model-based CLT reconstruction, radiative transfer equation (RTE) is used to mathematically describe the propagation of light in biological tissue. 40,41However, RTE has high computational complexity.To simplify the equation, di®usion approximation (DA) model is proposed and widely used. 16,42Based on ¯nite element analysis, a linear relation between the surface photon intensity and the distribution of the radioactive source is established, which can be de¯ned as follows: where È is the photon intensity of the surface.A is the system matrix determined by the characteristics of biological organs.X is the distribution of the radioactive source.

Deep learning-based CLT reconstruction
Deep learning-based CLT reconstruction uses neural network to establish the mapping relation between the photon intensity of surface and the radioactive source directly.Therefore, the deviation between the real photon propagation process and the simpli¯ed photon propagation process can be avoided.The deep learning-based method can be described as follows: where f nn denotes the reconstruction neural network.denotes network weight which is iteratively updated to reduce the error between the reconstructed radioactive source and the actual radioactive source.È denotes the surface photon intensity as an input to the network.X denotes the actual distribution of the radioactive source.

Three-dimensional graph convolution-based residual network for CLT
Generally speaking, the permissible domain can be discretized into n nodes, which can be generated by a FEM software.Each node has three-dimensional information as well as photon intensity information.However, three-dimensional information has been ignored in the previous deep learning-based methods.Hence, a residual network based on 3D graph convolution for CLT, named GCR-Net, which fully utilizes 3D information is described in detail in this section.

Three-dimensional graph convolution
As each node and its neighbors can be organized as graph structure, and inspired by DGCNN, 43 the EdgeConv layer is introduced into our network.In the EdgeConv layer, edge feature which describes the relationships between a node and its neighbors is proposed.It is de¯ned as follows: where e ij denotes edge feature.h denotes a nonlinear function with a set of learnable parameters.v i denotes the features of the node.v j denotes the features of its neighbors.To update the features of the node, an edge feature set of size k is calculated ¯rst.Then, an aggregation function such as mean operation is applied to aggregate edge features, which is de¯ned as follows: where v 0 i is the new features of the node.In Edge-Conv layer, the features of each node are updated.
To update the features, a graph structure is in need.In our network, each node was treated as a vertex in the graph.The input graph was constructed by executing a K-NN search to ¯nd the nearest neighbor in 3D coordinate space, and the k was set to 20.

Network architecture
GCR-Net is based on residual architecture.The basic residual unit in our model is named GCUnit.GCUnit consists of three EdgeConv layers.Each EdgeConv layer is followed by a leaky-ReLU layer.A skip connection is added between the input and the output of GCUnit.As shown in Fig. 1, GCR-Net is composed of ¯ve residual blocks, named GCBlock.There are two types of GCBlock.One consists of two GCUnits and the other consists of three GCUnits.A skip connection is introduced in GCBlock to reuse a previously learned feature, which connects from the input of GCBlock to its output.Before entering the ¯rst GCBlock, an EdgeConv layer with eight output channels is performed, followed by a leaky-ReLU layer.At the end of the last GCBlock, an EdgeConv layer with 1 output channel is performed and then a leaky-ReLU layer is attached.The network con¯guration we used is shown in Table 1.Taking GCBlock 1 as an example, it consists of two GCUnits.Each GCUnit consists of three EdgeConv layers with 8, 16, and 16 output channels, respectively.It should be noted that, in our network, normalized coordinates and intensity are taken as inputs, and the size is n Â 4. The output of GCR-Net is the reconstructed distribution of the radioactive source.

Implementation details and evaluation metrics
The training and test of GCR-Net were implemented using Pytorch and Python 3.7.All operation was performed on a personal computer with an AMD Ryzen 7 1700 Eight-Core Processor 3.00 GHz CPU and a NVIDIA GeForce GTX 1080 Ti GPU.The optimization function of GCR-Net was Adam optimizer with a learning rate of 0.001.Mean Square Error (MSE) Loss was adopted as a loss function.
The performance of GCR-Net was quantitatively evaluated with the following three metrics: Location error (LE), Dice coe±cient, and signal-to-noise ratio (SNR).
LE is de¯ned as the position error between the reconstructed radioactive source and the actual radioactive source.L r represents the barycenter of the reconstructed radioactive source and L 0 represents the barycenter of the actual radioactive source.It can be de¯ned as follows: Dice coe±cient is used to evaluate the similarity between the reconstructed source region and the actual source region.S rec represents the reconstructed source region and S ac represents the actual source region.It can be de¯ned as follows: The SNR is used to measure the visual quality.u 0 i represents the energy of the ith node in the reconstruction result, u i represents the energy of the ith node in the original numerical model.n is the total number of nodes in the numerical model.It can be de¯ned as follows:

Experiments and Results
In this section, the performance of GCR-Net was evaluated by performing numerical simulations and in vivo experiments.IVTCG and MFCNN methods were used as baselines because they represent one of the most commonly used reconstruction methods in  model-based methods and deep learning-based methods respectively. 35,37This section is structured as follows: The ¯rst part shows the reconstruction performance of our network in numerical simulation.It contains three experiments, including singlesource experiments, dual-source experiments, and experiments of testing the in°uence of the order of input.The second part presents the reconstruction ability of our method in in vivo experiments.

Data collection
A heterogeneous cylindrical phantom with a 10 mm diameter and a 30 mm height was designed to simulate the mouse body.The phantom consists of ¯ve kinds of organs: Heart, lungs, bone, liver, and muscle.The phantom is exhibited in Fig. 2. Optical parameters of all organs are presented in Table 2, which were obtained from the literature. 44he numerical phantom was discretized into a standard mesh with 4626 nodes by using COMSOL Multiphysics V5.6 (COMSOL Inc., USA) software, from which its coordinates can be obtained. 45The order of input was determined at the discretized process by COMSOL Multiphysics.The standard mesh is exhibited in Fig. 2(b).A small spherical radioactive source with 1 mm diameter was set in the model to present the tumor.The photon wavelength was set to 650 nm to ensure appropriate penetration.As a data-driven method, plenty of data is in need for training.As it is di±cult to collect a large amount of data from in vivo experiments, Monte Carlo method was used for data collection, because it can simulate the transmission process of photons.All the Monte Carlo simulation was conducted using Molecular Optical Simulation Environment (MOSE 2.3). 46We simulated a total of 221 single-source samples to make the radioactive source cover the permissible region as much as possible.3000 dual-source samples were assembled by randomly selecting and adding the corresponding data of two single-source samples, which were as follows: where È dual represents the surface photon of the assembled sample, È i represents the surface photon of the ith single-source sample.X dual and X i represent the distribution of the internal source of the assembled sample and the ith single-source sample, respectively.S n is the set of single-source samples.Overall, in this paper, 221 single-source samples and 3000 assembled dual-source samples have been simulated.80% of the dataset was randomly selected for training and the others were selected for testing.

Single-source simulation result
To quantitatively analyze the e®ect of our network for single-source simulation, the mean value and deviation of di®erent evaluation indexes were calculated, as exhibited in Table 3.As suggested from Table 3, GCR-Net obtained the lowest mean location error, highest mean Dice, and highest mean SNR.It demonstrated that the proposed network stably improved the accuracy in localization and  morphology recovery.Furthermore, the network obtained the minimum LE value of 0.11 mm and the maximum dice value of 0.79 in some cases, which further proved its advantages in terms of localization capability and morphology recovery ability.Two samples, named model 1 and model 2, were used to measure the single-source reconstruction result of di®erent methods.In model 1, the source was set at (À3, À4, 3) mm, while the source in model 2 was set at (À6, À5, 17) mm.The quantitative results are exhibited in Fig. 3, which show that our network still performs better in LE, Dice, and SNR. Figure 4 shows the 3D visualization as well as the 2D cross sections of reconstruction results obtained by GCR-Net, MFCNN, and IVTCG, respectively.
Two additional single-source experiments have been conducted to prove our network is insensitive to the order of input.Two samples with single source, named model3 and model4, were used.In this experiment, the order of input has been disarranged.MFCNN was performed as a comparative experiment under the same conditions.It is noteworthy that the previously trained network was directly used without further training.The quantitative results are exhibited in Table 4.It can be demonstrated from Table 4 that MFCNN is sensitive to the input order.If the order of input is disarranged, the MFCNN cannot work well.On the contrary, GCR-Net can successfully reconstruct the distribution of the radioactive source with high accuracy in all experiments.

Dual-source simulation result
Dual-source simulation experiments were carried out to further evaluate the performance of our methods on multi-source simulation.The quantitative results are exhibited in Table 5.The results in Table 5 revealed that our method still had its superiority in LE, Dice, and SNR, which presented a promising performance of reconstruction.However, compared with the single-source samples, the results of the dual-source samples were worse, which was probably caused by the raise in reconstruction complexity.
Two samples, named model 5 and model 6, were used to measure the dual-source reconstruction result of di®erent methods.In model 5, the centers of dual source were set at (À1, 3, 10) mm, (À3, À4, 22) mm, respectively.In model 6, the sources were set at (À2, À3.5, 4) mm, (À4.5, 0, 5) mm.The quantitative results are exhibited in Fig. 5, which illustrate a consistent trend with the former experiments.Figure 6 illustrates the visualization of the dual-source reconstruction results and suggests that both reconstructed radioactive sources have a clear shape and are easy to distinguish.
Two additional dual-source experiments have also been conducted to prove our network is insensitive to the order of input.Two samples with dual source, named model 7 and model 8, were used.The experimental setup was the same as the formal experiments on disarranging the order of input for single-source simulation.The quantitative results are shown in Table 6.It can be seen from Table 6 that in the dual-source experiment, when the order of input points changes, our network can achieve  high accuracy in the reconstruction, while MFCNN still cannot work well.

In vivo CLT reconstruction
The performance of our method was further evaluated with implanted small animal experimental data.The experimental data were from the CLT/ micro-CT dual-modal system on adult nude mouse.
The mouse was placed on the automatic rotating stage, which was rotating at 360 degrees with 1 interval to capture the images of CLI and X-ray projection.The CLI signal with wavelength of 650 nm and white light data were collected using a cooled high-sensitivity EMCCD (À80 , iXonEM Ultra 888, UK).The CT volume data was acquired by a micro-CT system (tube voltage of 40 kV p, tube current of 300 mA).The torso section of the mouse is shown in Fig. 7(a).We segmented the main organs of the mouse, including muscle, lung, heart, stomach, liver, and kidney, and integrated these organs into the mouse model.The optical parameters of each organ are consistent with those in the literature. 47A spherical glass vessel with a radius of 1 mm containing about 800 AE 50 Ci 18 F-FDG was used to mimic the lesion containing a radionuclide probe, which was implanted in the designed location.Animal care and protocols were approved by the Fourth Military Medical University Animal Studies Committee.The mouse was used twice for the single-source experiment and dual-source experiment, respectively.The position of the source was set at a di®erent place each time.
In single-source experiments, the source was implanted into the mouse at a speci¯c location (22, 20, 31.5)mm.In dual-source experiments, the sources were implanted at (20.5, 27.5, 16.5) mm, (14.5, 20, 36) mm, respectively.The CLI imaging obtained from each implanted model was mapped and registered to the surface of the mouse torso.GCR-Net, MFCNN, and IVTCG methods were used to carry out the single-source experiment and  dual-source experiment.The results are shown in Table 7 and Fig. 7, which prove the superiority of our network in CLT reconstruction for in vivo experiments.

Discussion and Conclusion
In this paper, a novel GCR-Net was proposed to achieve the accurate reconstruction of CLT.GCR-Net was built by stacking residual blocks based on 3D graph convolution.The main advantages of our network are summarized as follows: First, our network fully utilizes three-dimensional information and improves the performance of CLT reconstruction.Secondly, our network is insensitive to the order of input, which is not able to be satis¯ed by the traditional fully connected network.
To verify the performance of GCR-Net, singlesource and dual-source numerical simulations were implemented.IVTCG method and MFCNN method were adopted for comparison.Quantitative analysis demonstrated that the GCR-Net showed a better performance in source localization and morphology recovery.To further evaluate the practicability of GCR-Net, in vivo experiments were performed.Consistent with our numerical simulation experiments, GCR-Net obtained the most accurate reconstruction results of radioactive probe distribution in biological tissue.
It is worth noting that our network is insensitive to the order of input.Our method can achieve accurate reconstruction results even if the order is disarranged.It can be explained by the characteristic of the EdgeConv layer.The EdgeConv layer is designed to be invariant to the ordering of neighbors.In EdgeConv layer, the applied aggregation function is a symmetric function.Thus, the output of EdgeConv layer is invariant to the permutation of the input.
However, there are still some limitations of GCR-Net.When calculating the edge feature set, the memory cost will be signi¯cantly increased with a higher setting of k in KNN method.Besides, the utilization of standard mesh depends on certain software and the registration between standard mesh and actual structure probably brings extra error in CLT reconstruction.
In conclusion, GCR-Net is proposed for CLT reconstruction.It is based on the 3D graph convolution and has the advantage of being insensitive to the order of input, which can achieve accurate reconstruction in source localization and morphology recovery.We hope our network can assist CLT reconstruction e±ciently.

Con°icts of Interest
The authors declare that there are no con°icts of interest related to this article.

Fig. 1 .
Fig. 1.The network structure of GCR-Net.Coordinates and surface photon intensity are taken as input, and the output is the ¯nal reconstruction result.

Fig. 2 .
Fig. 2. The numerical phantom.(a) shows a single-source phantom, and (b) shows the standard mesh.

Fig. 7 .
Fig. 7. Comparison results for in vivo experiments.(a) presents the torso section and the result reconstructed by GCR-Net method, (b) and (c) are the results obtained by the MFCNN method and IVTCG method for the single-source experiment.(d)-(f) presents the results reconstructed by GCR-Net, MFCNN, and IVTCG methods for the dual-source experiment.The black circle indicates the true position and size of the source.

Table 2 .
Optical parameters of di®erent organs.

Table 3 .
Quantitative comparison of single-source simulation (mean AE SD).

Table 4 .
Results of the experiments with disarranged order of input for single source.

Table 6 .
Results of the experiments with disarranged order of input for dual-source experiment.