Muon Tomography imaging improvement using optimized scattering tracks data based on Maximum Likelihood Method
Xiao-Dong Wang, Kai-Xuan Ye, Yi Wang, Yu-Lei Li, Xie Wei, Ling-Yi Luo, Guo-Xiang Chen
aa r X i v : . [ phy s i c s . i n s - d e t ] J un Prepared for submission to JINST
XIV th Workshop on RESISTIVE PLATE CHAMBERS AND RELATED DETECTORS02.19-23,2018Puerto Vallarta, Mexico
Muon Tomography imaging improvement using optimizedscattering tracks data based on Maximum LikelihoodMethod
WANG Xiaodong, a YE Kaixuan, b , c WANG Yi, d LI Yulei, a , d WEI Xie, a Luo lingyi, a and CHENGuoxiang a a School of Nuclear Science and Technology, University of South China, Hengyang 421001, China b Institute of Plasma Physics Chinese Academy of Sciences, PO Box 1126, Hefei, Anhui 230031, China c University of Science and Technology of China, Hefei, Anhui 230026, China d Tsinghua University, Department of Engineering Physics, Beijing, 100084,China
E-mail: [email protected]
Abstract: P oint of colsest Approche algorithm (PoCA) based on the formalism of muon radiogra-phy using Multiple Coulomb scattering (MCS) as information source is previously used to obtainthe reconstruction image of high Z material. The low accuracy of reconstruction image is causedby two factors: the flux of natural muon and the assumption of single scattering in PoCA algo-rithm. In this paper, the maximum likelihood method based on the characteristics of Gaussian-likedistribution of muon tracks by MCS is used to predict the optimal track of outgoing muon. Thereceiver operating characteristic (ROC) and the localization ROC (LROC) are used as two analysismethods to evaluate the quality of reconstruction image. From the results of simulation, the perfectdiscrimination of longitudinal materials could be well achieved by maximum likelihood algorithmand the discriminate ratio that is predicted by the maximum likelihood method is about 15 % higherthan that of predicted by PoCA algorithm method. It is seen that the maximum likelihood methodcan greatly improve the accuracy of reconstruction image. Keywords:
Maximum Likelihood Scattering-Expectation Maximization, Muon tomography, GEANT4
ArXiv ePrint : 1234.56789(???) ontents
Cosmic-ray muon has high penetrability, which has been used for non-destructing imaging fordecades by imaging methods ranging from 2-D analysis using measurement of muon attenuationto 3D muon tomography based on Multiple Coulomb scattering (MCS) of muons as they acrossthe materials. The attraction of muon imaging technology is no manufactured radiation source,no artificial dose, and high sensitivity to high Z material ,such as special nuclear material (SNM)refs. [1]. The scattering information of muon can be recorded by a pair of position sensitivedetectors, such as Drift-Tube detector refs. [1, 2], Gas Electron Multiplier (GEM) refs. [3, 4], Multi-gap Resistive Plate Chamber (MRPC) refs. [5–7], which are placed on both sides of the detectedobject. The Point of Closet Approach (PoCA) firstly proposed by Los Alamos National Laboratory(LANL), is one of the most widely used reconstruction algorithms ,which has been investigated fordetection of high-Z material refs. [8–10] in last decades. The low accuracy of reconstruction imageis caused by two factors: the low flux of natural muon and the assumption of single scattering inPoCA algorithm.In this paper, the superiority and ability of the maximum likelihood scattering-expectationmaximization (MLS-EM) refs. [11, 12] algorithm,is presented. The characteristics of Gaussian-like distribution of muon tracks by MCS is used to predict the optimal track of outgoing muon.Maximum likelihood is widely used in reconstruction of image, however, a different iterativeExpectation Maximization (EM) algorithm is introduced to find the ML estimate of scatteringdensity profiles of material, which is efficient and flexible refs. [14] in medical image reconstruction.The images of three models of detected material by comparing the qulity of images with the result– 1 –f closet approach (POCA) is studied. And the quality of image is researched by the ROC (ReceiverOperation Characteristic) curve as well as the localization ROC (LROC) using the MATLABsoftware.
The principle of tomography is illustrated in Fig.1, where tomography refers to the reconstructionimage of object from projections taken from many different directions. M rays sample the objectcharacteristic function when they pass through the imaging area along the line. When the i th raypasses the imaging area, its sampling (or signal) can be observed. The relationship between a ray’ssampling and the discrete object characteristic function can be described as following expression: P i = Õ j w ij f j , (2.1)where, the weight w ij is the path length of the i th ray through the j th pixel or voxel. The reconstructedcharacteristic function ˆ f can be estimated by solving the system of linear equations (2.1),. Sincemuon tomography is based on the traditional method, there are still several assumptions to bemodified: • The cosmic ray muon has a broad angular distribution around zenith with limited flux. • The ray signal, the angle of MCS is, the multiple Coulomb scattering is stochastic with azero-mean Gaussian, and the actual distribution even has heavier tail. • the rough location could be realized, and the structure of high Z material is estimated roughly.. f f · · · f N − f N · · · f ( x, y ) y x P P P M − P M P i P i − w ij Figure 1 . (color online) The principle of traditional tomography. A discrete model of the object charac-teristic function is adopted by assuming uniform values within each pixel or voxel, denoted by the values f , f , · · · , f N . And the rays sampling can be observed. – 2 –s is illustrated in Fig.2, the scattering angle is exaggerated, the observed data D i of muon isthe scattering angle ∆ θ i . The x and y planes for each of the scattering angle ∆ θ i to add reconstructedinformation is used in muon tomography, D x , i ( D y , i ) as: D x , i = ∆ θ x , i = ( θ x , out − θ x , in ) i D y , i = ∆ θ y , i = ( θ y , out − θ y , in ) i . (2.2) x HΔx ΔθIncoming muonOutgoing muon
Figure 2 . (Color online)The trajectory of muon through the material , it can be described by the scatteringangle ∆ θ and displacement ∆ x . Here, the path length of ray is approximately equal to the thickness ofmaterial H . The conditional probability distribution of the observed data D i may be approximated asGaussian distribution[9] with a zero mean, which is defined as: P ( D i | λ ) = √ π | Σ i | / exp − D i Σ i ! . (2.3)Where Σ i is the variance; λ is the scattering density ditribution.The variance Σ i can be expressed as Σ i = p r , i Õ j L ij λ j , (2.4)where L ij is similar to w ij , which is the path length of the i th ray through the j th voxel, and p r , i isthe momentum ratio which is inversely proportional to i th muon momentum p i . When considering– 3 –he muon detector noise, the variance is refined as Σ i = C i + p r , i Õ j L ij λ j , (2.5)where, C i is the contribution of the detector noise. The expression of (2.1), (2.4) and (2.5) have asimilar form, while (2.4) and (2.5) could express the relationship between the variance of scatteringangle and the scattering density for the stochastic signal. In this simulation, the Monte Carlo simulation package GEANT4 is used to simulate the MTspectrometer based on position sensitive detectors. The sensitive area has a volume of × × m filled with work gas. To achieve the reconstruction result accordant with practical circumstances,the number of simulated muon is . × corresponding to 5 min of exposure in experiment andthe size of material cubes are × × cm .Fig.3 shows that a group of position sensitive detectors are located above the object, and onebelow. The position, angle and momentum of muon will be recorded in these detectors. Figure 3 . (Color online) the principle of muon tomography. The muon has three interactions a, b and c:transmission, energy loss and MCS. The points A,B,C and D are the positions of muon in these detectorsrespectively, and other information is shown in the figure 3.
The physical process of muon, shown in Fig.3, is described in GEANT4 of the physics listQGSP-BERT-HP.The spectrum of generated muon is established by the spectrum of cosmic raymuons in the range of 3 to 100 GeV, which obeys the empirical formula refs. [14]: dIdE × dcos θ = . E − . © « + . E cos θ GeV + . + . E cos θ GeV ª®®¬ , (2.6)where, θ is plane angle from vertical and E is the muon energy. Fig.4 shows that the probability density distribution of scattering angle ∆ θ can be approximatedas Gaussian distribution (as the red fitting curve shown) when muons are passing through the– 4 –luminum, iron, lead and uranium. The variance of scattering angle is the function of atomic Znumber. The results show that the RMS (mrad) are about 4.33 for aluminum, 10.04 for iron, 18.69for lead, and 25.24 for uranium, respectively.The PoCA is a simple geometric algorithm under the assumption of single scattering in eachevent. This algorithm can quickly discriminate the location and structure of object, especially forhigh Z material, which has been validated through experimental studies refs. [3, 8–10, 15]. −20 0 2000.020.040.060.080.1 Theta (milliradian) in Aluminum P r o b a b ili t y d e n s i t y −50 0 5000.010.020.030.040.05 Theta (milliradian) in Iron P r o b a b ili t y d e n s i t y −100 0 10000.0050.010.0150.020.025 Theta (milliradian) in Lead P r o b a b ili t y d e n s i t y −100 0 10000.0050.010.0150.02 Theta (milliradian) in Uranium P r o b a b ili t y d e n s i t y σ(Al) = 4.33 mrad σ(Fe) = 10.04 mradσ(Pb) = 18.69 mrad σ(U) = 25.24 mrad Figure 4 . (color online) The probability density distribution of scattering angle may be approximated as azero-mean Gaussian distribution. The RMS scattering are about 4.33 (mrad) for Al (upper-left), 10.04 for Fe(upper-right), 18.69 for Pb (lower-left), and 25.24 for U (lower-right). The blue histogram and red curve arethe simulated and fitting results, respectively.
Firstly, the scattering angle is in the order of milliradian. The track of each muon can becomputed as a straight line connecting the incoming and outgoing points. Secondly, it should beassumed that the scatter of each event only occurred once time on the closest point to the incomingand outgoing tracks, which is called as the PoCA point. But in the view of the three dimensionsspatial, the incident and scattered tracks may not be coplanar and not intersect at a point. The pointcloset to the each pair line (incoming and outgoing) are computed by solving a linear algebraicformulation and the midpoint of common perpendicular are taken as the PoCA point. Fig.5 is theschematic of PoCA algorithm in the voxellation of imaging area.Thirdly, the ray signal can be marked as s i and defined as: s i = q ( ∆ θ x , i + ∆ θ y , i ) = " ( θ x , out − θ x , in ) + ( θ y , out − θ y , in ) i . (3.1)Moreover, since muons come from different directions so that path length is different for eachmuon, we modify (2.4) to compute the scattering density estimate ˆ λ for reconstruction. λ j = p r σ θ L = Õ i : L i j , s ij p r , i M j L ij , (3.2)– 5 – igure 5 . (color online) The PoCA algorithm schematic. The blue and green lines are incident and exitingtracks respectively. The yellow cubes are voxels through which the muons pass, and the red points areassumed to be PoCA points. where, the ray signal s ij is: s ij = ( s i Poca v oxel alon g path except Poca v oxel (3.3) The new Maximun Likelihood Scattering-Expectation Maximization (MLS-EM) is designed furtherby using information of scattering angles, which distribute the scattering location along the ray trackinstead of assigning to the PoCA point according to probability statistics. In order to compute theray path lengths through voxels, the entry points to PoCA points to exit points are connected toestimate.The total function of muon data may be written as: P ( D | λ ) = Ö i P ( D i | λ ) , (3.4)where the conditional probability distribution P ( D i | λ ) is given by (2.3).The MLS estimation of the scattering density ˆ λ can be solved by maximizing the log likelihoodfunction: ˆ λ M LS = ar g max λ > λ Air LP ( λ | D ) = ar g max λ > λ Air Õ i − lo g ( Σ i ) − D i Σ i ! ! , (3.5)an EM algorithm to maximize the log likelihood function are developed to calculating the scatteringdensity estimate, which is more efficient and flexible than traditional method.The following MLS-EM update equation for the scattering density at the n th iteration arederived as: ˆ λ n + j , M LS − E M = median i : L i j , B nij . (3.6)– 6 –ince the measurements in x and y are independent, B ij can be computes as the average of B x , ij and B y , ij : B nx , ij ( B ny , ij ) = λ nj + D x , i ( D y , i ) L ij Σ i − L ij Σ i ! × p r , i ( λ nj ) . (3.7)More detailed information can be found in refs. [9]. Fig.6 shows the perspective view of three different scenes for comparison between PoCA and MLS-EM algorithms: the horizontal, diagonal and vertical objects. The number of incident muons is2 . × , corresponding to about 5 mins of exposure in the experiment. The voxels are the size of50 × × mm placed in the imaging area whose size is 2 × × m . (a) (b) (c) Figure 6 . (Color online) The perspective view of the simulated horizontal (a), diagonal (b) and vertical (c)scenes.
Fig.7 shows the comparison results of horizontal scene between PoCA and MLS-EM. In thereconstructed image, the voxels with scattering density in a range of [-0.5,5) are colored by green(the low-Z), [5,30) blue (the medium-Z), and [30, + inf . ] red (the high-Z), respectively. In theMLS-EM reconstruction, the maximum times of iteration is set to 100, and the initial value of voxelscattering density is set to be of air: λ = − refs. [6]. Furthermore, the averages of scatteringdensity in the voxels where objects appear are used to demonstrate the results in brackets. Fig.7(a,c)are the 2D and 3D reconstructed images using PoCA algorithm, respectively, which shows that thescattering density estimate of (Al, Fe, Pb, U) are (1.72, 7.75, 54.63, 57.83 mrad / cm ). Fig.7(b,d)are the 2D and 3D reconstructed images using MLS-EM algorithm, respectively, which shows thatthe estimation of scattering density are 1.27, 7.98, 52.68, 62.70 for Al, Fe, Pb, U,respectively.Fig.8 displays the comparison results of diagonal scene between PoCA and MLS-EM, re-spectively, which show that the scattering density of (Al, Fe, U) are (1.36, 7.45, 63.61) by PoCAreconstruction and (1.73, 7.38, 68.37)by MLS-EM correspondingly. Fig.9 shows the comparisonresults of vertical scene between PoCA and MLS-EM, respectively, which reflect that the scatteringdensity of (Al, Fe, U) are (1.66, 9.96, 58.59) by PoCA reconstruction and (1.53, 8.11, 61.41)– 7 – (mm) Y ( mm ) −300 −200 −100 0 100 200 300−300−200−1000100200300 20406080100120X(mm) Y ( mm ) −300 −200 −100 0 100 200 300−300−200−1000100200300 20406080100120−500−400−300−200−1000 100 200 300 400 500−500−400−300−200−1000100200300400500−500−400−300−200−1000100200300400500 X(mm)Y(mm) Z ( mm ) −500−400−300−200−1000 100 200 300 400 500−500−400−300−200−1000100200300400500−500−400−300−200−1000100200300400500 X(mm)Y(mm) Z ( mm ) (a) (b)(c) (d) Al(1.72) Fe(7.75)Pb(54.63) U(57.83) Al(1.27) Fe(7.98)Pb(52.68)U(62.70)
Figure 7 . (color online) The comparison reconstruction in 2D (a,b) and 3D (c,d) of the horizontal materialbetween PoCA (a,c) algorithm with MLS-EM (b,d) algorithm. by MLS-EM correspondingly.Compared with the PoCA algorithm, the MLS-EM algorithm canreconstruct better location and appearance, particularly in vertical scenes.
X(mm) Y ( mm ) −300 −200 −100 0 100 200 300−300−200−1000100200300 20406080100120 −500−400−300−200−1000 100 200 300 400 500−500−400−300−200−1000100200300400500−500−400−300−200−1000100200300400500 X(mm)Y(mm) Z ( mm ) (a)(c) X(mm) Y ( mm ) −300 −200 −100 0 100 200 300−300−200−1000100200300 20406080100120−500−400−300−200−1000 100200 300400 500−500−400−300−200−1000100200300400500−500−400−300−200−1000100200300400500 X(mm)Y(mm) Z ( mm ) (d)(b) Al(1.36)U(63.61) Fe(7.45) Fe(7.38)U(68.37) Al(1.73)
Figure 8 . (color online) The comparison reconstruction in 2D (a,b) and 3D (c,d) of the diagonal materialbetween PoCA (a,c) algorithm with MLS-EM (b,d) algorithm.
The goal of muon tomography is to discriminate and exclude the presence of high-Z material whichcan achieve a high efficiency and keeping false positives rate low in a short reconstructed time. Weplot the ROC curve and the localization ROC (LROC) curve that has been commonly used in thebinary discrimination system to evaluate the image quality in PoCA and MLS-EM reconstruction.Three sets of samples are generated, as are shown in Fig.10. A set of sample hiding the target (atungsten cube inside the iron volume) is regarded as the "signal" sample, the other two sets of samplewithout the target (empty or iron cube inside the volume) are regarded as the "noise" samples. The50 images per set are tested when the event of muon is 8 . × or less. The maximum of scattering– 8 – Z ( mm ) (c) X(mm) Y ( mm ) −300 −200 −100 0 100 200 300−300−200−1000100200300 20406080100120140−500−400−300−200−1000 100 200300 400500−500−400−300−200−1000100200300400500−500−400−300−200−1000100200300400500 X(mm)Y(mm) Z ( mm ) (b)(d) X(mm) Y ( mm ) −300 −200 −100 0 100 200 300−300−200−1000100200300 20406080100120140160 (a) Al(1.66)Fe(9.96)U(58.59) Al(1.53)Fe(8.11)U(61.41)
Figure 9 . (color online) The comparison reconstruction in 2D (a,b) and 3D (c,d) of the vertical materialbetween PoCA (a,c) algorithm with MLS-EM (b,d) algorithm.
Figure 10 . (Color online) The ROC (LROC) analysis illustration for muon tomography. density , λ max , in the image is compared with a selected threshold, T . The true positive rate (TPR,sensitivity) is considered as the probability to trigger the alert in the "signal" image when λ max > T .Otherwise, the false positive rate (FPR, 1-specificity) is defined as the probability to exclude thepresence of target in the "noise" image when λ max < T . A series of pair of sensitive and specificitycan be obtained by varying the threshold. The sensitivity plotted against 1-specificity is the ROCcurve, the perfect method would yield a point in the upper-left corner of the ROC curve. In thatcase, the area under the curve (AUC) is equal to 1.To consider the localization, the LROC curves are also performed, which have been used inthe medical imaging community for assessing the lesion localization performance refs. [12]. Thetest statistic in LROC differ in that the true positive appears as if and only if the λ max is obtainedat the target location in the "signal" sample. Fig.11 exhibits the analysis of ROC and LROC curvebetween PoCA and MLS-EM algorithm. It is clear that PFR, MLS-EM can achieve much higherTPR than PoCA algorithm. For discriminating the signal sample from the noise sample, the MLS-EM reconstruction achieve AUC of 0.9964 in the air noise, 0.9764 in the iron noise for the ROCcurve, 0.858 in the air noise, and 0.8592 in the iron noise for the LROC curve, respectively. Theresults are much higher than those of PoCA reconstruction (0.8708, 0.8366, 0.7312, and 0.798– 9 –orrespondingly). The comparison reflects that the MLS-EM algorithm can significantly improvethe performance of ROC and LROC to the PoCA approach. S e n s iti v it y ( T P R ) MLS−EM,AUC=0.9964PoCA,AUC=0.8708 S e n s iti v it y ( T P R ) MLS−EM,AUC=0.9764PoCA,AUC=0.8366 S e n s iti v it y ( T P R ) MLS−EM,AUC=0.858PoCA,AUC=0.7312 S e n s i t i v i t y ( T P R ) MLS−EM,AUC=0.8592PoCA,AUC=0.798 (a) (b)(c) (d)
Figure 11 . (Color online) The analysis of ROC (LROC) and AUC between PoCA and MLS-EM reconstruc-tion: (a) ROC: W vs. Air, (b) ROC: W vs. Fe, (c) LROC: W vs. Air, (d) LROC: W vs. Fe.
This paper describes two tomographic reconstructions based on MCS of natural muons. ThePoCA algorithm is employed to detect the high-Z material. And considering the actual application,a instance of EM method is designed to solve large size ML problems which differs and moreefficient. The simulated results suggest that, MLS-EM algorithm can significantly improve theimaging performance of muon tomography.
Acknowledgments
This work was supported by National Natural Science Foundation of China (Grant No.11605086),the foundation of Natural Science Foundation of Hunan Province (Grant No. 2018JJ3422), the thefoundation of Education Department of Hunan Province (Grant No. 15B205), and the foundationof open project of state key laboratory of particle detection and electronics.)
References [1] Borozdin K N, Hogan G E, Morris C, et al,
Surveillance: Radiographic imaging with cosmic-raymuons , Nature . (2003) 277.[2] Bittner B, Dubbert J, Horvat S, et al, Development of fast high-resolution muon drift-tube detectorsfor high counting rates, Nucl. Instr. and Meth. A. (2011) 154-157. – 10 –
3] K. Gnanvo, P. Ford, J. Helsby et al,
Performance expectations for a tomography system using cosmicray muons and micro pattern gas detectors for the detection of nuclear contraband , Proceeding ofIEEE Nuclear Science Symposium Conference Record, Dresden, Germany . (2008) 1278-1284.[4] Jared Sturdy,
Technical Report CMS-CR-2016-047 , CERN, Geneva . (2016)[5] Ye J, Cheng J, Yue Q, et al,
Studies on RPC position resolution with different surface resistivity ofhigh voltage provider , Nuclear Science Symposium Conference Record, 2008, NSS’08, IEEE . (2008)917-918.[6] Baesso P, Cussans D, Thomay C, et al,
Toward a RPC-based muon tomography system for cargocontainers , Journal of Instrumentation . (2014) C10041.[7] Wang X, Zeng M, Zeng Z, et al, The cosmic ray muon tomography facility based on large scaleMRPC detectors , Nucl. Instr. Meth. A , (2015) 390-393.[8] Miyadera H, Borozdin K N, Greene S J, et al, Imaging Fukushima Daiichi reactors with muons , AipAdvances. (2013) 052133.[9] J. O. Perry, Advanced applications of cosmic-ray muon radiography , PhD thesis (The University ofNew Mexico) . (2013)[10] Rhodes C J,
Muon tomography: looking inside dangerous places, Science progress. (2015)291-299.[11] AP Dempster, NM Laird, and DB Rubin, Maximum likelihood from incomplete data via the EMalgorithm , PJR Stat. Soc. B. (1977) 1-38.[12] Green P J, Latuszynski K, Pereyra M, et al, Bayesian computation: a summary of the current state,and samples backwards and forwards , Statistics and Computing. (2015) 835-862.[13] Schultz L J, Blanpied G S, Borozdin K N, et al, Statistical reconstruction for cosmic ray muontomography, IEEE transactions on Image Processing. (2007) 1985-1993.[14] Nakamura K, Particle Data Group, Review of particle physics, Journal of Physics G: Nuclear andParticle Physics. (2010) 075021.[15] Schultz L J, Cosmic ray muon radiography , PhD thesis (Portland State University) . (2003). (2003)