A fuzzy artificial neural network-based method for Cerenkov luminescence tomography

Cerenkov Luminescence Tomography (CLT) is a non-invasive three-dimensional in vivo detection technology. However, due to the illposedness of CLT, the reconstructed result has many artifacts, which will mislead the researchers to make a wrong diagnostic decision. Enlightened by the development of artificial neural networks, we proposed a Fuzzy Autoencoder Clustering method to eliminate these artifacts and improve the reconstruction quality. To assess the performance of our method, several numerical simulation experiments and real physical phantom experiments are conducted. Compared with the raw reconstruction results and the commonly used manual threshold processed ones, it is demonstrated that our method is capable of filtering the artifact areas effectively, making reconstruction results clearer. It is anticipated that the method presented in this paper will help advance the CLT technology and promote the clinic translation of CLT technology. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5088234


I. INTRODUCTION
Cerenkov luminescence imaging (CLI) is a promising optical molecular imaging (OMI) tool due to its advantages of relatively short acquisition time and wide clinical radionuclide probes available. [1][2][3][4] The basic principle of CLI is that when charged particles (such as β+ positrons produced by radionuclides) travel faster than light in the medium, spectrally continuous radiation (also named as Cerenkov Luminescence, CL) will be emitted, which is known as the Cerenkov effect. 5 Since its first application to the detection of the distribution of radionuclide in the small animal by Roberson et al. in 2009, CLI has attracted an increasing number of researchers' attention and become a hot research topic in the field of OMI. 6 Compared to other imaging modalities in OMI filed, the biggest advantage of CLI is its clinical-used radionuclide probes approved by the Food and Drug Administration (FDA). 7,8 However, CLI can only obtain two-dimensional distribution of radionuclide probe and lacks depth information, 9 so in 2010, a tomographic imaging technology based on CLI image was proposed by Li et al. and named as Cerenkov Luminescence Tomography (CLT). 10 Nevertheless, for CLT two major problems also exist. First, it is known that the reconstruction problem of CLT is highly illposed, which means that it is difficult to accurately reconstruct the distribution of radionuclides. This will produce some artifacts in the recovered image. [11][12][13] Second, the CL emitted from the tissue surface is weak in real experiments, which leads to the fact that the noise signal may easily affect the collected CL data. This defect will make the reconstruction result easily interfered by noise, which limits the application of CLT technology. [14][15][16] To address these problems, researchers have proposed various methods to improve the reconstruction accuracy, such as Tikhonov regularization, sparsity regularization, ln regularization method, multispectral data-based method, multi-modal imaging system, etc. [17][18][19] However, from the reconstruction result, it can be found that not only the nodes of the light source have a value larger than 0, but also many other nodes have such values due to the ill-posedness of CLT, forming artifacts which will mislead researchers to make a wrong diagnostic decision. The most frequently-used method in CLT study is to adopt a hard threshold to eliminate these artifacts. The threshold is typically set as 70% of the reconstructed maximum energy value, and then the nodes with the value lower than that are set to 0. 20 But in some cases, these artifacts have a similar value to that of the light source, which makes it difficult to set an appropriate threshold. Therefore, an effective postprocessing method is needed for CLT. In 2018, Yi et al. proposed an adaptive threshold method for Fluorescence Molecular Tomography (FMT) technology. 21 In the method, the reconstructed nodes are sorted according to value and then a threshold is calculated. However, the situation that some nodes in artifacts have a relatively large value (that may be similar to that of light sources) is still not considered.
In this paper, we proposes a new reconstruction method for CLT based on artificial neural network (ANN). First, the value of a node and its neighbors in the reconstruction result are assembled by Near Four Element Fields (NFEF) strategy. Then, the assembled data are used to train a Fuzzy Autoencoder (FAE) network. In this step, fuzzy number is introduced to enhance the robustness and the capability of the network, consequently reducing the impact of noise on the results. Finally, the latent features obtained by FAE are used as the inputs of the Fuzzy C-means (FCM) clustering algorithm to cluster the reconstructed nodes into two groups, and the group with the lower total energy will be set to 0. 22 This paper is mainly divided into four sections. Section II gives the corresponding theoretical knowledge of the method, as well as the details of all the steps. Section III gives out the experimental settings and the corresponding results as well as their analysis. The conclusion and discussion are presented in Section IV.

A. CLT reconstruction
For CLT, Diffusion Approximation (DA) model can be used to describe the transmission of photons. The DA model combined with the Robin boundary condition is described in Eq. (1): where Φ(r) denotes the luminescent flux density of the photon at the point r in the region Ω, D(r) = 1 3 µa + µ ′ s is the diffusion coefficient, µa and µ ′ s represent the absorption coefficient and the reduced scattering coefficient, respectively, ξ represents the point at the boundary of the surface of the object, Z is determined by the degree of matching of the corresponding biological tissue and the refractive index of the air, and n is the unit normal vector whose direction points to the outside of ∂Ω.
With the finite element method (FEM) framework, a linear relationship between the photon energy flow rate of the outer surface and the unknown internal source distribution can be established as Eq. (2): 19 where Am×n is a system matrix describing the transmission of Cerenkov radiation in biological tissue, Φ is the surface photon flow rate, and Xn ×1 is the unknown source distribution to be solved. In this paper, the Incomplete Variables Truncated Conjugate Gradient (IVTCG) algorithm 23 is employed as the reconstruction algorithm.
In the CLT reconstruction methods available, the manual threshold is often used to filter the results. Empirically, the threshold is set as 70% of the maximum energy value Ψ shown as Eq. (3).
This approach is simple in practice, but in some cases, it may engender biased results because the mere artificial threshold setting is lacking in flexibility, leading to erroneous conclusions. Hence, focusing on this issue, we introduce the FAE-FCM into post processing.
B. FAE-FCM method for post-processing 3) Clustering: applying FCM clustering algorithm to Y to obtain the effective energy Ee and the background noise En respectively, and the final result vector X ′ being determined based on Ee.

C. Data processing
In order to apply X to train the AE network, and inspired by the definition of Markov Random Field (MRF), 24 we proposed an NFEF strategy to increase the numbers of each row of X. In NFEF strategy, the energy value of each node and the energy of its four nearest where the first column of the M matrix is the original X vector, xNi,j represents the j nearest neighbors of the i th element in the original X.
To improve the robustness of the training process of AE, the fuzzy number, which is flexible is introduced, and thus the FAE is got. 25 Given a real number k and an input i, there will be a mapping relationship: for i = k, f (i) = 1. When the absolute value of i approaches to infinity, f (i) = 0. Here, f (i) is the corresponding membership function. In this paper, the symmetric triangular fuzzy number will be adopted: , where β is a positive constant that controls the width of ⌢ k. Actually, ⌢ k corresponds to a closed interval [kL, kR] and its center is k.
For a given constant σ ∈ [0, 1], we can define a new membership function: The new mapping relation is called the σ-cut of is also a fuzzy number which has the same center as ⌢ k does. We can also define its left bound In the latter case, the fuzzy number degrades to a common real number. Thus, we can fuzzy each element of M matrix in Eq. (4) to get: The loss function of AE is used to measure the proximity of the input data and the output. However, a fuzzy number corresponds to an interval that contains infinitive points and we cannot compute the reconstructed error for every point. So it is necessary to discretize the continuous fuzzy domain to facilitate the establishment of loss function, and the σ-cut is used here. In this paper, the closed interval [0, 1] is divided into m equivalent subintervals and thus we have m + 1 endpoints, and then we define Where m = 10 in this paper. According to the previous description, it can be seen that each σ ∈ [0, 1) corresponds to two real numbers M[σ] L and M[σ] R , and Eq. (6) can be then expressed as: In this paper, Eq. (8) is split to the left and the right as the input of AE, whose rationality has been demonstrated by our previous work: 26

D. Feature extraction
After the data processing, the original result X has been fuzzified, and two fuzzy matrices in Eq. (9) are obtained as the training data of the AE network. The AE consists of an input layer, an output layer and a hidden layer, is used to find suitable parameters ξ 1 and ξ 2 to satisfy the relationships between input SL, SR and output SL ′ , SR ′ as: 27,28 Moreover, the AE can be seen the combination of an encoder and a decoder, which can be seen in Eq. (11) and Eq. (12), separately.
is the Sigmoid activation function. It should be noted that in this paper, the left and right parts of AE network are both 40 hidden neurons. And in the implementation of our network, according to the previous work, W ′ can be directly replaced by W T . 29 Then, the loss function L in this paper is established as: After minimizing the loss function, one parameter set, can be determined and the output of encoder: y 1 and y 2 are used as the latent features of the input data. 26 The basic AE network is shown in Fig. 2(a), while the FAE used in this paper is depicted in Fig. 2(b).

E. Clustering
After the feature extraction process, we get the latent features of the input. Then, we combine y 1 and y 2 to generate matrix Y, as shown in Fig. 1. Then, FCM is selected as the clustering algorithm. 22 Given U = (µi,j)n×c is a fuzzy classification matrix, where c is the number to be classified (in this paper, c=2), and µi,j is the membership of sample j for classification, and the sample type Yj is the feature of the hidden layer of the FAE. Construct the objective function as follows: 22 where m ∈ (1, ∞) is the fuzzy index, which determines the degree of classification matrix ambiguity. In this paper, m = 2, and di,j = ||Yj − ωi||, ωi denotes the center of the cluster i.
Using the Lagrangian multiplier method to minimize objective function, the following membership optimization expression can be obtained: And the optimized cluster center: After clustering, feature matrix Y can be clustered into two new matrices. As each row of the matrix Y corresponding to one element in vector X, thus, the two new matrices represent the latent features of the two parts of X. We define the part with higher sum value is the effective nodes, and the corresponding matrix is named as Ee, while the other one named as En. So, after setting the value of the nodes corresponding to En to 0, we can get the final reconstruction result vector X ′ .

III. EXPERIMENT AND RESULTS
In this paper, the reconstruction method is implemented by using MATLAB 2018A. The main configuration of the computer is a 3.70 GHz Intel Core i7 CPU, 16 GB RAM, and NVIDIA GTX 1080Ti GPU.
In order to evaluate the performance of the proposed method in this paper, several sets of numerical simulations and physical phantom experiments are conducted. Several quantitative evaluation indicators including the centroid error (EC), the weighted center error (EWC), the reconstructed volume ratio (RV ), and the coincident tetrahedral ratio (R T ) are employed to quantitatively assess the reconstruction results.
The centroid error EC is defined as the Euclidean distance between the centroid coordinates of the reconstruction result (x, y, z) and the true Cerenkov source (x 0 , y 0 , z 0 ).
The weighted center error EWC is defined as the Euclidean distance between the weighted center coordinates of the reconstruction result (xc, yc, zc) and the center coordinates of the real Cerenkov source (xc 0 , yc 0 , zc 0 ).
The reconstruction volume ratio RV is defined as the ratio between the total volume of the reconstructed tetrahedron Vrec and the volume of the real source Vtrue, as follows: The reconstruction algorithm obtains the energy values of a group of nodes, which will form the corresponding internal tetrahedral structure. Accordingly, these tetrahedral structures constitute the final reconstruction area. The overall volume of the tetrahedron can be calculated by the coordinates of the four nodes of the tetrahedron. But, some nodes of some tetrahedrons may have no energy value, which means that only a part of the volume of the tetrahedron belongs to the reconstruction result. For these tetrahedrons, the volume corresponding to each tetrahedron can be equally divided into four parts, and if there are several nodes whose energy values are non-zero, it is considered that the tetrahedron contributes the volume of the corresponding number of parts in the reconstruction area. For the mesh used in this paper, this method can effectively calculate the total volume of the reconstruction tetrahedron Vrec.
The coincident tetrahedron ratio R T is the number of intersections between the index of the reconstructed tetrahedron (Irec) and the index of real source tetrahedron (Iture), divided by the number of  unions between the index of the reconstructed tetrahedron and the index of real source tetrahedron, as follows:

A. Numerical mouse simulation
In the following simulations, the numerical mouse model is used. 30 The optical parameters are given in Table I. 18 Spheres with a radius of 0.8mm serve as the light sources. In the single source simulation, one light source is placed at (12mm, 14mm, 25mm), as shown in Fig. 3(a). In the double source simulation, two light sources are placed at (9mm, 7mm, 18.5mm) and (26mm, 7mm, 18.5mm), separately, as shown in Fig. 3(b). The Monte Carlo (MC) method is used to generate the photon flow rate distribution of the forward model. And the FEM method is used to discretize the numerical mouse into 2206 surface triangular patches, 47407 internal tetrahedrons, and 8588 mesh nodes.

ARTICLE
scitation.org/journal/adv  Table II. According to the experimental results, it can be seen that although the artifacts with lower energy values can be filtered by the manual threshold, those with higher values still exist. Compared with the manual threshold, the proposed FAE-FCM method yields better results, which is due to the combination of the energy distribution of the light source and the implicit spatial geometry information of nodes.  Table III further presents the quantitative analysis of the results. When the nodes with higher energy value are not located inside the light source, the artificial filtering method will yield a deteriorated result. Fig. 5 shows the results obtained by different methods intuitively, and a comparison of Figs. 5(b) and (c) will show the advantage of the proposed FAE-FCM framework over the artificial threshold approach. The quantitative results in Table III also show the error evaluating indicators (EC and EWC) of the manual threshold method increase significantly.

Double source numerical mouse simulation
To further test the proposed method with more complicated settings, two more experiments of double source numerical mouse simulation are conducted. In the first experiment, there are two sources with different power density in the numerical mouse model (one is 1 nW/mm 3 and the other is 0.5 nW/mm 3 ), and the location of light sources are (9mm, 7mm, 18.5mm) and (26mm, 7mm, 18.5mm). Fig. 6(a) is the forward simulation result, while Fig. 6(b) is the reconstruction result. It can be concluded that the two light sources can be well reconstructed. For the second experiment, the initial values of light

ARTICLE
scitation.org/journal/adv sources are the same, but the distance between them is closer than that in the previous experiment, with the location being (9mm, 7mm, 18.5mm) and (13mm, 8mm, 18.5mm). Fig. 6(c) is the forward simulation result, while Fig. 6(d) is the reconstruction result.

Materials and instruments
Here, 18 F-deoxyglucose ( 18 F-FDG) is selected as the radioactive material. It is made by GE's MINItrace cyclotron. The radiochemical purity of 18 F-FDG is 97%, radioactive concentration is 5 GBq/mL, pH is 8.0, and the purity is greater than 99.9%.  6. (a) and (c) is the forward simulation result, while (b) and (d) is the reconstruction result. In (a), the power densities of the two sources are different; and in (c), the distance between the two sources is closer than that of (a), with the densities are the same.
The optical information of the simulated surface was collected by IVIS imaging system manufactured by Calipers Life Sciences, USA.
In this section, a cube and a cylinder phantom made of polyoxymethylene are used. The height of the cube and the cylinder are 25mm, and the length of the cube and the diameter of the cylinder are both 25mm. A cylindrical hole with a diameter of 3mm and depth of 12.5mm is drilled at the top of two phantoms respectively. By injecting 10µg of 18F-FDG into the hole, a cylindrical light source about 2mm high with a center coordinate of (6.25mm, 0mm, 1mm) was constructed. After injecting, the top of the phantom holes is sealed with masking tape to avoid interference from other factors. The aperture of phantom was sealed with a light-shielding tape. The phantoms geometry and the source setting are shown in Figs Table IV. It can be seen from Fig. 8 that the all of reconstructed source size is significantly reduced compared to the numerical simulation, which mainly caused by the inevitable loss and error in signal acquisition. This leads to an over-filtering phenomenon in Figs. 8(b) and (e). A large amount of the lower energy portion is filtered out, resulting in a too small reconstructed result. Similar conclusions can be drawn from

Result of cylinder phantom experiment
The comparison results among the raw reconstruction result of the cylinder phantom experiment, the processed result with the manual threshold and the result with the FAE-FCM method are shown in Fig. 9. The quantitative analysis of the results is shown in Table V. Similar to cubic phantom experiment, the higher energy portion of the reconstructed result is located outside the 18 F-FDG source, as shown in Figs. 9(a) and (d). This leads to that most of the energy information inside the source is filtered out with the    Table V. On the contrary, the proposed post-processing method performed better. It filters out the portion located in the non-source region even though the mesh nodes in this portion have a relatively higher energy value.

IV. CONCLUSION AND DISCUSSION
In this paper, we propose a novel post-processing method for CLT reconstruction. Compared to the traditional threshold method, the proposed method is free from artificial threshold, and performs better even when the nodes with the maximum intensity value in the reconstructed source region are not located inside the true ones. The reconstruction results of simulations and physical phantom experiments presented in Section III demonstrate that the FAE-FCM post-processing method can improve the quality of the reconstructed images by reducing artifacts and increasing the location accuracy.
However, it should be noted that in the algorithm of this paper, there are two parameters that may have a significant impact on the reconstruction result, and they are the unit number of hidden layer neurons and the maximum fuzzy factor. In all of the experiments involved in this paper, the left and right parts of AE network are both 40 hidden neurons. We use 80 as the unit number of hidden layer neurons because a larger number cannot improve the performance of the reconstruction result, according to an additional experiment shown in Fig. 10(a), which depicts the reconstruction results obtained by using a serious of the different number of hidden layer neurons, and other parameters are consistent with the numerical experiments of single light source. It is obvious that when the number of hidden layer neurons is larger than 80, the result will not change significantly. As for the selection of maximum fuzzy factor (10 is used in this paper), we have also conducted a serious of experiments to investigate its impact on the reconstruction result. Four experiments with different maximum fuzzy factors are conducted and the results are shown in Fig. 10(b), and other parameters are consistent with the cubic physical phantom experiment, we can see that when the maximum fuzzy factor exceeds 10, the results mainly remain unchanged.
As one main contribution of this work is that the FAE-FCM is a post-processing part of the CLT framework, it can be combined with any other reconstruction method. So, we hope this work may provide a useful tool for the CLT, as well as other optical molecular tomography technologies.