Hybrid model for muon tomography and quantitative analysis of image quality

Hybrid model for muon tomography and quantitative analysis of image quality 1 first-author https://orcid.org/0000-0002-2402-2066 1 1 https://orcid.org/0000-0003-4006-5282 1 1 1 1 https://orcid.org/0000-0002-3165-6964 1 1 corresp *wangxd@usc.edu.cn *wangxd@usc.edu.cn wangxd@usc.edu.cn https://orcid.org/0000-0001-9838-1195 School of Nuclear Science and Technology, University of South China, Hengyang 421001, China School of Nuclear Science and Technology, University of South China, Hengyang 421001, China Muon tomography is a novel method for the non-destructive imaging of materials based on muon rays, which are highly penetrating in natural background radiation. Currently, the most commonly used imaging methods include muon radiography and muon tomography. A previously studied method known as coinciding muon trajectory density tomography, which utilizes muonic secondary particles, is proposed to image low and medium atomic number ( Z ) materials. However, scattering tomography is mostly used to image high- Z materials, and coinciding muon trajectory density tomography exhibits a hollow phenomenon in the imaging results owing to the self-absorption effect. To address the shortcomings of the individual imaging methods, hybrid model tomography combining scattering tomography and coinciding muon trajectory density tomography is proposed and verified. In addition, the peak signal-to-noise ratio was introduced to quantitatively analyze the image quality. Different imaging models were simulated using the Geant4 toolkit to confirm the advantages of this innovative method. The simulation results showed that hybrid model tomography can image centimeter-scale materials with low, medium, and high Z simultaneously. For high- Z materials with similar atomic numbers, this method can clearly distinguish those with apparent differences in density. According to the peak signal-to-noise ratio of the analysis, the reconstructed image quality of the new method was significantly higher than that of the individual imaging methods. This study provides a reliable approach to the compatibility of scattering tomography and coinciding muon trajectory density tomography. Monte Carlo simulation Muon tomography Image reconstruction Introduction 1 Cosmic ray muons are high-energy charged particles produced by the cascade clustering of primary cosmic rays with atoms in the atmosphere as they are emitted toward the Earth's surface [ 1 1 ] . With an average energy of 3–4 GeV and a lifetime of approximately 2.2 μs, their mass is approximately 207 times heavier than that of an electron, and their flux at sea level is approximately 1 cm -2 ·min -1 [ 2 2 ] . The advantages of muon tomography include high penetration, the existence of muons everywhere, and its non-destructive effects on humans, which allow it to be widely applied in various cross-disciplinary studies. Currently, there are two main types of cosmic ray muon imaging methods: transmission imaging and scattering tomography. Transmission imaging technology reconstructs images based on the change in muon intensity before and after the muon penetrates a material. Based on this, transmission imaging mainly focuses on objects with large volumes. In 1970, Alvarez applied this method to the field of archaeology and discovered dark chambers inside the pyramids [ 3 3 ] . This method has subsequently been widely used for active volcanoes, underground tunneling [ 4 4 , 5 5 ] , hydrogeological studies [ 6 6 ] , and exploration of unknown chambers within the Pyramids of Khufu [ 7 7 ] . Muon scattering tomography (MST) employs the multiple Coulomb scattering properties of cosmic muons in a material, which reconstructs the structure and outline of a high- Z material by changing the angle of incident trajectories and outcoming trajectories. Los Alamos National Laboratory (LANL) first introduced this method in 2003 and applied it to identify high- Z materials in a short time [ 8 8 ] . LANL successively raised several imaging algorithms, such as the point of closest approach (PoCA), maximum likelihood scattering and displacement (MLSD), and maximum a posteriori (MAP). [ 9 9 , 11 11 ] In the experimental concept, LANL also proposed the first prototype of muon tomography based on a drift tube structure, which is able to confirm whether there is nuclear material hidden in the detection area [ 12 12 ] . However, owing to its small sensitive area, the imaging results of this prototype had clear defects. Therefore, different types of position-sensitive detectors were designed to image different materials [ 13 13 - 17 17 13-17 ] . A wide range of MST applications have been demonstrated, ranging from environmental protection to homeland security. This method was subsequently applied to nuclear security studies, such as reactor monitoring [ 18 18 , 19 19 ] and container detection [ 20 20 ] . An essential application of muon scattering tomography was performed by a research team in the United States to reconstruct images of nuclear material littered in a reactor after the Fukushima nuclear power plant accident [ 21 21 ] . Compared with transmission imaging, scattering tomography possesses better imaging accuracy and less imaging time, but it is usually used for high- Z and high-density materials and is not capable of imaging low -Z materials. Both scattering tomography and transmission imaging inherently infer the information of the target material through the change in muon trajectories or flux before and after they interact with materials. In fact, the secondary particles generated from the interaction between the muon and the material also carry significant information about the target object. If muonic secondary particles are fully applied, more abundant information about the object will be obtained. There are two main methods of reconstructing objects using muonic secondary particles. The first is to mark the incident muon using the secondary muonic neutron produced in special nuclear materials, which can confirm the existence of special nuclear materials in a highly penetrating and low-dose manner without disclosing its factual information [ 22 22 - 24 24 22-24 ] . However, because the number of muons captured to produce secondary neutrons is small, the imaging time of this method is too long, and the imaging quality is relatively poor. The second method is based on the muon and muonic secondary particle coincidence detection method proposed by Mrdja [ 25 25 - 27 27 25-27 ] , which was used to experimentally reconstruct a three-dimensional image of a bovine femur, achieving the first muon imaging of low- Z matter. The research team at the University of South China also performed related studies and theoretically proposed a novel four-dimensional (4D) imaging method using the coincidence detection technique of cosmic-ray muons and their secondary particles [ 28 28 ] . This novel 4D imaging method was verified in a Monte Carlo simulation, demonstrating the function of secondary electrons and gamma rays in the process of tomography and broadening the application field of muon tomography. However, this method requires a longer imaging time and mainly aims at imaging low- Z materials, and its reconstruction quality for high- Z materials is not promising. To obtain the structure and location information of high-, medium-, and low- Z materials in a single experiment, a muon imaging method that is applicable to complex situations and compatible with multiple types of information was studied. To address the shortcomings of individual imaging methods, we propose a muon hybrid model tomography that combines scattering tomography and coinciding muon trajectory density tomography to achieve high-quality three-dimensional (3D) reconstructed images and acquire an excellent ability to distinguish between high-, medium-, and low- Z materials. The peak signal-to-noise ratio (PSNR) [ 29 29 ] value was then introduced to quantitatively analyze hybrid model tomography. The results of the analysis showed that the quality of the reconstructed image was clearly better than that of scattering tomography and trajectory density tomography. Simulation parameters 1 Geometry of the detector system 2 A three-section removable hybrid model tomography system (combination of sections (a), (b), and (c)) based on the Geant4 toolkit [ 30 30 ] was constructed, as shown in Fig. 1 Fig. 1 , which was composed of a trajectory detection module (combination of sections (a) and (c)) and secondary particle signal trigger module (combination of sections (a) and (b)). Through this construction, scattering tomography and trajectory density tomography could be combined. Sections (a) and (c) are the three layers of the muon trajectory detectors set at the top and bottom of the detection area, respectively. The function of the trajectory detector is to record the position of the muon hitting the detectors to fit the incident and outgoing tracks. Section (b) consists of four plastic scintillators with a box-like shape that can generate a trigger signal when secondary particles are incident. The trajectory detector measured 500 mm × 500 mm × 10 mm, and the interval between the two detectors was 50 mm. As the existing spatial resolution of the Micromegas detector can reach ~65 μm [ 31 31 ] , the spatial resolution of the detector was set to 100 μm, which can fulfill the demand for accuracy in the experiment. For the secondary particle signal trigger module, the individual detector was described as a plastic scintillator detector of C 9 H 10 according to the EJ200 model, with dimensions of 500 mm × 500 mm × 50 mm. The distance between each section was 100 mm. Previous results acquired by our research team [ 28 28 ] indicated that the energy of the secondary particles is sufficiently large, and part of the energy deposited in the plastic scintillator is higher than the baseline of the detectors and other nuclear electronic plugs. As a result, it is reasonable to distinguish secondary particles from background noise. Therefore, the detection efficiency was assumed to be 100% during the simulation. During each event, the value of the muon trajectory detection module and secondary particle signal trigger module was set to “0” until secondary particles passed through the scintillation detector, then this value was reset to “1”. Simulation package and muon generator 2 Included in the Geant4 toolkit was a particle source package to generate an incident muon, a system modeling module to construct a detector and imaging object, a particle physics process list, and a data output port, which can simulate the movement of muons and their secondary particles in the object and position-sensitive detectors. Subsequently, the muon location and direction were acquired to fit the muon trajectory. The incident particle source was assigned to the independent physics program Cosmic-ray Shower Library (CRY) [ 32 32 ] . Concept of imaging methods 1 Concept of coinciding muon trajectory density tomography 2 Coinciding muon trajectory 3 Muons are high-energy particles that lose energy mostly through ionization, bremsstrahlung, pair electron production, and photonuclear interactions. According to a study by Bogdanov et al. [ 33 33 ] , muons mainly interact with materials through ionization when their energy is less than 10 GeV. When a muon is incident, the trajectory detection module responds and records its position information. Subsequently, the muon knocks out an electron of an atom of the material and forms a δ electron. The δ electron will then generate photons through the bremsstrahlung process, and a part of the photons will produce electrons via the three main ways in which photons interact with materials. From the viewpoint of low- or medium- Z materials, these secondary particles have a probability of overcoming the self-absorbing effect of the material and escaping from its surface. At different values of Z , the ratio of the number of escaped secondary photons to the number of escaped secondary electrons is different. The escaped secondary particles enter the scintillation detector to produce a signal. Therefore, simulation results indicate that muonic secondary particles primarily consist of electrons and gamma rays. A method for imaging target objects based on coinciding muon trajectories is proposed in this paper. The principle of this imaging method is as follows: In an ideal situation, if the trajectory detection module and secondary particle signal trigger module respond simultaneously within a time gate of approximately 1–2 ns, the incident muon trajectory will be classified as a coinciding muon trajectory because it has a high likelihood of passing through the target object. The thickness of the object along a certain direction can be inferred using a trajectory density reconstruction algorithm when a large number of coinciding muon trajectories accumulate in that direction. If sufficient trajectories are obtained, a 3D image of the object can be reconstructed. However, there are two experimental situations that prevent the accurate acquisition of all secondary particles. The first is when an incident muon produces several secondary particles and hits the same scintillator detector, making it difficult to distinguish several secondary particles owing to the time resolution of the scintillation detector. The second is when several secondary particles hit various detectors and generate a recordable signal. Because they are inside a time gate, they are all experimentally treated as one signal. Consequently, if an incident muon generates many secondary particles, the number of muon trajectories will only increase by one throughout the screening process. Therefore, optimization to increase the number of effective counts is proposed, the method of which involves recording the signal once there is one plastic scintillator that coincides with the trajectory-detection module. In addition, we assumed that the coinciding muon trajectory is a straight line. Imaging algorithm 3 An algorithm based on this principle is illustrated in Fig. 2 Fig. 2 . First, as shown in Fig. 2(a) Fig. 2(a) , the imaging area is divided into a series of mathematical planes, and every coordinate of every cross point between all coinciding muon trajectories and all planes is calculated. When the mathematical planes are perpendicular to the Z -axis, the method of calculating the coordinates of the cross points is equivalent to the method used to confirm the intersection of a 3D line with a plane that is perpendicular to the Z -axis. Subsequently, the area with the highest cross-point density is chosen, as shown in Fig. 2(b) Fig. 2(b) . The cross points outside the target object are represented by yellow dots, and the red dots represent the cross points inside the object, which can be considered an effective reconstruction event. Finally, the object is divided into small voxels, and the chosen points in the voxels are accumulated to assign the value of the voxels, as shown in Fig. 2(c) Fig. 2(c) . By displaying the value of every voxel, a two-dimensional (2D) image of the object is reconstructed. In 3D imaging, it is important to confirm the height of the object. To obtain altitude information, we assume that the volume surrounded by four plastic scintillators is a cube and set the interval of its top and bottom surfaces as the altitude detection area. Subsequently, a checkout of the distribution of cross points along the Z -axis is performed to determine where a sudden increase or decrease in distribution occurs. The distribution of cross points suddenly changes when the planes are higher than the topmost section of the object or lower than its bottommost section. In this way, the height of the object is obtained through the distance between the two abruptions, and a 3D reconstruction is performed. Concept of scattering tomography 2 Multiple Coulomb scattering occurs when a muon penetrates the material. During this process, the muon deflects several times at small angles, and as the deflection angle accumulates, large angular and positional deflections of the outgoing trajectory compared with the incident trajectory occur. The total deviation angle has an approximately Gaussian distribution. Eq. (1) provides the standard deviation of the deflection angle, and Eq. (2) yields the radiation length X 0 , which is a property of the material. <math display="block" id="M1"><mrow><msub><mi>σ</mi><mi>θ</mi></msub><mo>≈</mo><mrow><mo>[</mo><mrow><mn>13.6</mn><mo>/</mo><mrow><mo>(</mo><mrow><mi>β</mi><mo>⋅</mo><mi>C</mi><mo>⋅</mo><mi>p</mi></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>⋅</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>L</mi><mo>/</mo><msub><mi>X</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow><mrow><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup><mo>⋅</mo><mrow><mo>[</mo><mrow><mn>1</mn><mo>+</mo><mn>0.038</mn><mrow><mo>(</mo><mrow><mi>L</mi><mo>/</mo><msub><mi>X</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M2"><mrow><msub><mi>X</mi><mn>0</mn></msub><mo>=</mo><mn>716.44</mn><mi>A</mi><mo>/</mo><mrow><mo>[</mo><mrow><mi>ρ</mi><mo>⋅</mo><mi>Z</mi><mrow><mo>(</mo><mrow><mi>Z</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>⋅</mo><mi>ln</mi><mrow><mo>(</mo><mrow><mn>287</mn><mo>/</mo><msup><mi>Z</mi><mrow><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> where p is the muon's momentum, L is the thickness passed by the muon in the target matter, β is the muon's velocity in relation to the speed of light, C is the speed of light, and A , ρ , and Z are the atomic mass number, density, and atomic number of the target material, respectively. The distribution of the scattering angles is linked to the density and atomic number of the material, as shown in Eqs. (1) and (2). Assuming that the muon is monoenergetic, the scattering angle after Coulomb scattering is proportional to the atomic number of the material in muon imaging detection; the higher the atomic number, the larger the Coulomb scattering angle. Therefore, the distribution of the scattering angles is capable of distinguishing high-, medium-, and low-Z materials. Concept of hybrid model tomography 2 To compensate for the defects of individual imaging methods and improve the quality of the reconstructed images, a combination of trajectory density tomography and scattering tomography is proposed, which is known as muon hybrid model tomography. First, two images are obtained through scattering and trajectory density tomography. The 2D image is then reconstructed based on the coefficient distribution image hybrid technique, as shown in Fig. 3 Fig. 3 . For 3D images, the target object is divided into several layers based on the computed tomography (CT) principle. Each layer is treated with reference to the 2D reconstruction method, after which the processed data are imported into the 3D Slicer toolkit. The coefficient distribution image hybrid technique assigns a coefficient according to the grayscale value and the detailed information of the original image. This technique has the advantages of fast calculation speed, ability to increase the signal-to-noise ratio, and ability to maintain the details of the original images. The discrete 2D wavelet transform method is introduced to achieve this coefficient distribution image hybrid technique. The discrete 2D wavelet decomposition and reconstruction process for individual images are shown in Fig. 3 Fig. 3 . The decomposition process can be described as follows: First, discrete 1D wavelet decomposition is performed on every column of the original image to acquire the low-frequency component L and high-frequency component H , as given by <math display="block" id="M3"><mrow><mi>L</mi><mrow><mo>[</mo><mrow><mi>m</mi><mo>,</mo><mi>n</mi></mrow><mo>]</mo></mrow><mo>=</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>K</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>K</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mi>x</mi></mstyle><mrow><mo>[</mo><mrow><mi>m</mi><mo>,</mo><mn>2</mn><mi>n</mi><mo>−</mo><mi>k</mi></mrow><mo>]</mo></mrow><mo>⋅</mo><mi>g</mi><mrow><mo>[</mo><mi>k</mi><mo>]</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M4"><mrow><mi>H</mi><mrow><mo>[</mo><mrow><mi>m</mi><mo>,</mo><mi>n</mi></mrow><mo>]</mo></mrow><mo>=</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>K</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>K</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mi>x</mi></mstyle><mrow><mo>[</mo><mrow><mi>m</mi><mo>,</mo><mn>2</mn><mi>n</mi><mo>−</mo><mi>k</mi></mrow><mo>]</mo></mrow><mo>⋅</mo><mi>h</mi><mrow><mo>[</mo><mi>k</mi><mo>]</mo></mrow><mo>,</mo></mrow></math> where m and n are the height and width of the original image, respectively. In addition, g ［ k ］ and h ［ k ］ represent the low-pass and high-pass filters, respectively. The main component of an image is low-frequency information, which comprises the basic grayscale of the image. High-frequency information forms the edges and details of the image. Discrete 1D wavelet decomposition is then performed for every row of the acquired low-frequency component L to obtain the low-frequency (LL) and high-frequency (LH) components along the horizontal direction. Every row of the acquired high-frequency component H is also subjected to this decomposition to obtain the low-frequency (HL) and high-frequency (HH) components along the horizontal direction. The formulas used to calculate LL, LH, HL, and HH are similar to Eqs. (3) and (4). After completing the above steps on the original image, two groups of LL, LH, HL, and HH (transformed data) are obtained to merge the images reconstructed from scattering tomography and trajectory density tomography. Then, if the discrete 2D wavelet reconstruction is performed, the fused image will be achieved. The discrete 2D wavelet reconstruction process for the individual images is reproduced below. By performing an inverse discrete 1D wavelet transform on every column and every row of the transformed data, the reconstructed image is obtained. Discrete wavelet decomposition is a process that decomposes the signal according to the low frequency and directional high frequency, which can decompose the wavelet to meet the demand for image accuracy. Results for simple model 1 The results of the three imaging methods are presented in Fig. 4 Fig. 4 . The target objects were uranium, iron, and water cubes, all with a volume of 50 mm × 50 mm × 50 mm from left to right. The number of incident muons was 7 × 10 6 , which was equal to approximately two days of actual imaging time according to the natural flux of the muons. A 3D image reconstructed using trajectory density tomography is shown in Fig. 4(a) Fig. 4(a) . According to the results, it is clear that there was a void inside the imaging results of uranium due to the strong self-absorption effect of uranium, which caused a lower emission of secondary particles from the inside than from the edge. In contrast, the self-absorbing effect of water is weaker than that of other materials, but its cross-section for producing secondary particles is relatively smaller, which also caused a slight hollow phenomenon. However, the water cube was still visible, and a more accurate image could be obtained by increasing the number of incident muons. The reconstructed image of iron had the best imaging quality and contained no voids. Therefore, the method was most suitable for imaging low- Z and medium- Z materials, with the best results for medium Z materials, and the reconstructed images of high- Z materials had voids in the middle. The reconstructed results of scattering tomography shown in Fig. 4(b) Fig. 4(b) indicate that scattering tomography can better obtain the shape and location of uranium and acquire a high-quality image. Moreover, the method can be used to obtain the approximate position of the iron block. However, the completely reconstructed outline image was not clearly discriminated, which could be improved by increasing the imaging time. In addition, this method is incapable of inspecting images of water. Consequently, the scattering signal of the high- Z material was stronger than that of the low- and medium- Z materials, and the lower the Z , the weaker the signal. Therefore, scattering tomography is primarily sensitive to high- Z materials. However, its disadvantage of being unable to reconstruct low- Z materials was also apparent. A comparison with the image reconstruction results obtained using the individual imaging methods revealed that hybrid model tomography can reconstruct images of uranium, iron, and water simultaneously (see Fig. 4(c) Fig. 4(c) ). The shape and location information of the three items were acquired with clear differentiation. When imaging the uranium block, the voids also disappeared, achieving a more accurate image reconstruction. It can be observed from the results that hybrid model tomography can solve the hollow phenomenon when imaging high- Z materials and compensate for the fact that scattering tomography cannot reconstruct low- Z materials. Reconstructed images with better quality can clearly identify high-, medium-, and low- Z materials. Results for complex models 1 To further verify the advantages of hybrid model tomography, two complex models containing high-, medium-, and low- Z materials were constructed: a Rubik's cube model and cylindrical model (see Fig. 5 Fig. 5 ). The cylindrical model used to study the influence of shielding on different imaging techniques was composed of uranium, iron, and water with diameters of 40 mm, 80 mm, and 120 mm, respectively, and a height of 120 mm. The Rubik's cube model used to study the discriminating ability of materials with similar Z was composed of nine tiny cubes, each with a volume of 30 mm × 30 mm × 30 mm and a spacing of 5 mm. This model contained more materials to simulate objects in an actual situation. The uranium, copper, and calcium oxide were arranged in the first row. Aluminum, iron, and lead were placed in the second row, and water, iron, and tungsten were placed in the third row. The iron in the third row worked as a reference group for the second row, which was used to study the influence of position on the self-absorbing effect. The reconstruction area was 1200 mm × 1200 mm with 5 mm × 5 mm pixels. The imaging results for the cylindrical model are shown in Figs. 6(a) and (b) Figs. 6(a) and (b) . The results of scattering tomography demonstrated that the method can accurately image both iron and uranium materials, and uranium wrapped in the center can be visualized from the top view. However, the water surrounding the surface could not be imaged. As shown in the imaging results, trajectory density tomography can reconstruct water, iron, and uranium simultaneously but cannot distinguish the uranium wrapped in the center. This is because the secondary particles generated by the muon interaction with the material were absorbed in the process of penetrating the shielding material, resulting in a reduction in the coinciding muon trajectories of the uranium material. Figures 6(c) and (d) Figures 6(c) and (d) show the reconstructed images of the Rubik's cube model. The results revealed that scattering tomography is mainly sensitive to high- Z materials and can roughly distinguish Pb and U of different atomic numbers in high- Z materials. Nevertheless, the imaging results for medium- Z materials were poor, and the imaging of low- Z materials was impracticable. Therefore, it was difficult to reconstruct the complete Rubik's cube model. Trajectory density tomography can determine the general shape of the cubes in the Rubik's cube model and distinguish between 5 mm gaps. However, the imaging results for high- Z materials still exhibited hollow phenomena, and the self-absorption effect was more pronounced for the iron block in the center than for the reference iron block. The imaging results of hybrid model tomography for both the models are shown in Figs. 6(e) and (f) Figs. 6(e) and (f) . It can image the three materials of uranium, iron, and water in the cylindrical model concurrently and provide superior 3D imaging results, as well as identify shielded uranium. For the Rubik's cube model, hybrid model tomography can clearly distinguish materials with close atomic numbers whose densities differ significantly, such as Pb and U, Pb and W in high atomic numbers. It can better reconstruct the cubes of high- and medium- Z materials and obtain location information on low- Z materials. In addition, it can roughly distinguish the 5 mm interval between cubes with high and medium Z in the model. Thus, hybrid model tomography can compensate for the disadvantages of individual methods, increase the image quality, improve the discriminating ability for high atomic number materials with close Z values, and precisely image high-, medium-, and low- Z materials simultaneously, which guarantees its potential in the imaging of complex structural objects. Quantitative analysis of image quality 1 To quantitatively evaluate the reconstructed image quality of hybrid model tomography, this section introduces PSNR analysis. The PSNR is the most common and widespread standard for objectively evaluating the quality of images and works on the basis of the error between the pixels of different images, which is to be evaluated by Eq. (5). In other words, the PSNR is an error-sensitive image quality evaluation method. The mean square error (MSE) represents the mean square error between the reconstructed image and the reference image, which is given by Eq. (6). <math display="block" id="M5"><mrow><mi>P</mi><mi>S</mi><mi>N</mi><mi>R</mi><mo>=</mo><mn>10</mn><mo>⋅</mo><msub><mrow><mi>log</mi></mrow><mrow><mn>10</mn></mrow></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mi>M</mi><mi>A</mi><msubsup><mi>X</mi><mi>I</mi><mn>2</mn></msubsup></mrow><mrow><mi>M</mi><mi>S</mi><mi>E</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>=</mo><mn>20</mn><mo>⋅</mo><msub><mrow><mi>log</mi></mrow><mrow><mn>10</mn></mrow></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mi>M</mi><mi>A</mi><msub><mi>X</mi><mi>I</mi></msub></mrow><mrow><msqrt><mrow><mi>M</mi><mi>S</mi><mi>E</mi></mrow></msqrt></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M6"><mrow><mi>M</mi><mi>S</mi><mi>E</mi><mo>=</mo><mfrac><mn>1</mn><mrow><mi>m</mi><mi>n</mi></mrow></mfrac><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mrow><msup><mrow><mrow><mo>‖</mo><mrow><mi>I</mi><mrow><mo>(</mo><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow><mo>)</mo></mrow><mo>−</mo><mi>K</mi><mrow><mo>(</mo><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msup></mrow></mstyle></mrow></mstyle><mo>,</mo></mrow></math> where m and n represent the height and width of the image to be analyzed, respectively. I ( i, j ) denotes the pixels in the reference image, and K ( i, j ) indicates the pixels in the image to be evaluated. MAX I represents the maximum color value of the two images. The reference image and the image to be evaluated in this study are all gray images stored in the form of nonlinear 8-bit sampling pixels; that is, there are 256 layers of grayscale in the images. The PSNR is measured in decibels (dB), and the larger the value, the less distortion between the assessed image and the reference image. The “UI” letter model (see Fig. 7(a) Fig. 7(a) ) was formed of cubes with 30 mm side lengths sewed together, whose material was set as lead. The letter “U” was composed of seven cubes, and the letter “I” was composed of three cubes. The pixel size was 5 mm × 5 mm, and different numbers of incident muons were converted into various imaging times. The reconstruction results of the hybrid model, scattering, and trajectory density tomography at various imaging times are shown in Fig. 7(a) Fig. 7(a) . Within one hour, the reconstructed image distortion was severe, making it impossible to acquire shape data for the “UI” letter. As the imaging time increased to 4 h, trajectory density tomography still could not reconstruct this model, whereas scattering and hybrid model tomography could roughly reconstruct its outline; however, the image quality of hybrid model tomography was better. Although all three imaging methods gathered information on the object's position and shape over time, hybrid model tomography was clearly superior to the other imaging methods, as shown visually. The PSNR values obtained by comparing the reconstructed images of the three imaging methods with the reference image for various imaging times are shown in Fig. 7(b) Fig. 7(b) . Hybrid model tomography had the highest PSNR value at the same imaging time, quantitatively indicating that its reconstructed image was better than that of the individual methods. In addition, a new method to quantitatively evaluate whether there is a split on the surface of an object is proposed. Taking the “UI” letter model as an example, we studied the changing trend of the gray value for each pixel along one certain direction (red line shown in Fig. 7(c) Fig. 7(c) ) and used it as a criterion for judging. If there is indeed a split, the gray value of that area will decrease heavily, and the closer the gray value is to “0”, the better the judging effect. The results revealed that hybrid model tomography can collect the letter model's location and shape information in 4 h and can discern a 10 mm gap. Conclusion 1 Muon hybrid model tomography was proposed by organically combining muon scattering tomography and trajectory density tomography. The imaging ability of these methods for materials with different atomic numbers was investigated by setting different imaging models. As shown in the imaging results, scattering tomography was more suitable for reconstructing the images of high- Z materials. However, the reconstruction results of this method for materials with a lower Z were poor, particularly for low- Z materials. For trajectory density tomography, the imaging effect of medium Z and low Z was better than that of scattering tomography. However, for high- Z materials, such as uranium and lead, because of the self-absorbing effect, the reconstructed image exhibited a hollow phenomenon. According to the imaging results of hybrid model tomography, this method could distinguish different materials in addition to reconstructing high-, medium-, and low- Z materials on an hourly scale. It was capable of discriminating materials used in shielding and reconstructing shielding objects to avoid misjudgment caused by individual imaging methods. For materials with similar atomic numbers, hybrid model tomography clearly distinguished those with obvious differences in density. Subsequently, a PSNR analysis was introduced to evaluate the quality of the reconstructed image. Finally, a new method was proposed to quantitatively analyze whether there was a split in the object. It can be seen from the results of this study that the longer the imaging time, the better the quality of these imaging methods. The acquired image of hybrid model tomography was apparently better than that of the two individual imaging methods and could obtain the position and shape information of the imaging model as well as distinguish 10 mm intervals. References S.H. Neddermeyer , C.D. Anderson , Cosmic-ray particles of intermediate mass . Phys. Rev. 54 , 88 ( 1938 ). doi: 10.1103/PhysRev.54.88.2 http://doi.org/10.1103/PhysRev.54.88.2 Cosmic-ray particles of intermediate mass M. Tanabashi , P.D. Grp , K. Hagiwara et al. , Review of particle physics: Particle data group . Phys. Rev. D. 98 , 030001 ( 2018 ). doi: 10.1103/PhysRevD.98.030001 http://doi.org/10.1103/PhysRevD.98.030001 Review of particle physics: Particle data group L.W. Alvarez , J.A. Anderson , F.E. Bedwei et al. , Search for hidden chambers in the pyramids . Science. 167 , 3919 ( 1970 ). doi: 10.1126/science.167.3919.832 http://doi.org/10.1126/science.167.3919.832 Search for hidden chambers in the pyramids V. Tioukov , A. Alexandrov , C. Bozza et al. , First muography of Stromboli volcano . Sci. Rep. 9 , 6695 ( 2019 ). doi: 10.1038/s41598-019-43131-8 http://doi.org/10.1038/s41598-019-43131-8 First muography of Stromboli volcano R. Han , Q. Yu , Z. Li et al . Cosmic muon flux measurement and tunnel overburden structure imaging . J. Instrum. 15 , P06019 ( 2020 ). doi: 10.1088/1748-0221/15/06/P06019 http://doi.org/10.1088/1748-0221/15/06/P06019 Cosmic muon flux measurement and tunnel overburden structure imaging R. Nishiyama , A. Ariga , T. Ariga et al. , First measurement of ice-bedrock interface of alpine glaciers by cosmic muon radiography . Geophys. Res. Lett. 44 , 12 ( 2017 ). doi: 10.1002/2017GL073599 http://doi.org/10.1002/2017GL073599 First measurement of ice-bedrock interface of alpine glaciers by cosmic muon radiography K. Morishima , M. Kuno , A. Nishio et al. , Discovery of a big void in Khufu's Pyramid by observation of cosmic-ray muons . Nature. 552 , 387 - 401 ( 2017 ). doi: 10.1038/nature24647 http://doi.org/10.1038/nature24647 Discovery of a big void in Khufu's Pyramid by observation of cosmic-ray muons K.N. Borozdin , G.E. Hogan , C. Morris et al. , Surveillance: Radiographic imaging with cosmic-ray muons . Nature. 422 , 277 - 279 ( 2003 ). doi: 10.1038/422277a http://doi.org/10.1038/422277a Surveillance: Radiographic imaging with cosmic-ray muons S. Tripathy , J. Datta , A. Jash et al. , Material identification with cosmic ray muons using RPCs . J. Instrum. 14 , C07007 ( 2018 ). doi: 10.1088/1748-0221/14/07/C07007 http://doi.org/10.1088/1748-0221/14/07/C07007 Material identification with cosmic ray muons using RPCs S. Xiao , W.B. He , M.C. Lan et al . A modified multi-group model of angular and momentum distribution of cosmic ray muons for thickness measurement and material discrimination of slabs . Nucl. Sci. Tech. 29 , 28 ( 2018 ). doi: 10.1007/s41365-018-0363-7 http://doi.org/10.1007/s41365-018-0363-7 A modified multi-group model of angular and momentum distribution of cosmic ray muons for thickness measurement and material discrimination of slabs C.M. Liu , Q.G. Wen , Z.Y. Zhang et al . Study of muon tomographic imaging for high-Z material detection with a Micromegas-based tracking system . Radiat. Detect. Technol. Methods 4 , 263 - 268 ( 2020 ). doi: 10.1007/s41605-020-00179-9 http://doi.org/10.1007/s41605-020-00179-9 Study of muon tomographic imaging for high-Z material detection with a Micromegas-based tracking system C.L. Morris , C.C. Alexander , J.D. Bacon et al. , Tomographic imaging with cosmic ray muons . Sci. Glob. Secur. 16 , 37 - 53 ( 2008 ). doi: 10.1080/08929880802335758 http://doi.org/10.1080/08929880802335758 Tomographic imaging with cosmic ray muons X. Wang , M. Zeng , Z. Zeng et al. , The cosmic ray muon tomography facility based on large scale MRPC detectors . Nucl. Instrum. Meth. A 784 , 390 ( 2015 ). doi: 10.1016/j.nima.2015.01.024  http://doi.org/10.1016/j.nima.2015.01.024  The cosmic ray muon tomography facility based on large scale MRPC detectors Z.W. Pan , J.Y. Dong , X.J. Ni et al. , Conceptual design and update of the 128-channel μSR prototype spectrometer based on musrSim . Nucl. Sci. Tech. 30 , 123 ( 2019 ). doi: 10.1007/s41365-019-0648-5 http://doi.org/10.1007/s41365-019-0648-5 Conceptual design and update of the 128-channel μSR prototype spectrometer based on musrSim X.Y. Pan , Y.F. Zheng , Z. Zeng et al. , Experimental validation of material discrimination ability of muon scattering tomography at the TUMUTY facility . Nucl. Sci. Tech. 30 , 120 ( 2019 ). doi: 10.1007/s41365-019-0649-4 http://doi.org/10.1007/s41365-019-0649-4 Experimental validation of material discrimination ability of muon scattering tomography at the TUMUTY facility X.L. Qian , H.Y. Sun , C. Liu et al. , Simulation study on performance optimization of a prototype scintillation detector for the GRANDProto35 experiment . Nucl. Sci. Tech. 32 , 51 ( 2021 ). doi: 10.1007/s41365-021-00882-2 http://doi.org/10.1007/s41365-021-00882-2 Simulation study on performance optimization of a prototype scintillation detector for the GRANDProto35 experiment Z.Y. Li , Y.M. Zhang , G.F. Cao et al. , Event vertex and time reconstruction in large-volume liquid scintillator detectors . Nucl. Sci. Tech. 32 , 49 ( 2021 ). doi: 10.1007/s41365-021-00885-z http://doi.org/10.1007/s41365-021-00885-z Event vertex and time reconstruction in large-volume liquid scintillator detectors P.L. Rocca , V. Antonuccio , M. Bandieramonte et al. , Search for hidden high- Z materials inside containers with the Muon Portal Project . J. Instrum. 9 , C01056 ( 2013 ). doi: 10.1088/1748-0221/9/01/C01056 http://doi.org/10.1088/1748-0221/9/01/C01056 Search for hidden high-Z materials inside containers with the Muon Portal Project J. Armitage , J. Botte , K. Boudjemline et al. , First images from the cript muon tomography system . Int. J. Mod. Phys. Conf. Ser 27 , 1460129 ( 2014 ). doi: 10.1142/S201019451460129X http://doi.org/10.1142/S201019451460129X First images from the cript muon tomography system G. Blanpied , S. Kumar , D. Dorroh et al. , Material discrimination using scattering and stopping of cosmic ray muons and electrons: Differentiating heavier from lighter metals as well as low-atomic weight materials . Nucl. Instrum. Meth. A 784 , 352 - 358 ( 2015 ). doi: 10.1016/j.nima.2014.11.027 http://doi.org/10.1016/j.nima.2014.11.027 Material discrimination using scattering and stopping of cosmic ray muons and electrons: Differentiating heavier from lighter metals as well as low-atomic weight materials K. Borozdin , S. Greene , Z. Lukic et al. , Cosmic ray radiography of the damaged cores of the Fukushima reactors . Phys. Rev. Lett. 109 , 152501 ( 2012 ). doi: 10.1103/physrevlett.109.152501 http://doi.org/10.1103/physrevlett.109.152501 Cosmic ray radiography of the damaged cores of the Fukushima reactors E. Guardincerri , J. Bacon , K. Borozdin et al. , Detecting special nuclear material using muon-induced neutron emission . Nucl. Instru. and Meth. A 789 , 109 - 113 ( 2015 ). doi: 10.1016/j.nima.2015.03.070 http://doi.org/10.1016/j.nima.2015.03.070 Detecting special nuclear material using muon-induced neutron emission T.B. Blackwell , V.A. Kudryavtsev , Simulation study into the identification of nuclear materials in cargo containers using cosmic rays . J. Instrum. 10 , T04003 ( 2015 ) doi: 10.1088/1748-0221/10/04/T04003 http://doi.org/10.1088/1748-0221/10/04/T04003 Simulation study into the identification of nuclear materials in cargo containers using cosmic rays H.Y. Lan , T. Song , J.L. Zhang et al. , Rapid interrogation of special nuclear materials by combining scattering and transmission nuclear resonance fluorescence spectroscopy . Nucl. Sci. Tech. 32 , 84 ( 2021 ). doi: 10.1007/s41365-021-00914-x http://doi.org/10.1007/s41365-021-00914-x Rapid interrogation of special nuclear materials by combining scattering and transmission nuclear resonance fluorescence spectroscopy I. Bikit , D. Mrdja , K. Bikit et al. , Novel approach to imaging by cosmic-ray muons . Europhysics Lett. 113 , 58001 ( 2016 ). doi: 10.1209/0295-5075/113/58001 http://doi.org/10.1209/0295-5075/113/58001 Novel approach to imaging by cosmic-ray muons D. Mrdja , I. Bikit , K. Bikit et al. , First cosmic-ray images of bone and soft tissue . Europhysics Lett. 116 , 48003 ( 2016 ). doi: 10.1209/0295-5075/116/48003 http://doi.org/10.1209/0295-5075/116/48003 First cosmic-ray images of bone and soft tissue G. Galgoczi , D. Mrdja , I. Bikit et al. , Imaging by muons and their induced secondary particles—a novel technique . J. Instrum. 15 , C06014 ( 2020 ). doi: 10.1088/1748-0221/15/06/C06014 http://doi.org/10.1088/1748-0221/15/06/C06014 Imaging by muons and their induced secondary particles—a novel technique X.T. Ji , S.Y. Luo , Y.H. Huang et al. , A novel 4D resolution imaging method for low and medium atomic number objects at the centimeter scale by coincidence detection technique of cosmic-ray muon and its secondary particles . Nucl. Sci. Tech. 33 , 2 ( 2022 ). doi: 10.1007/s41365-022-00989-0 http://doi.org/10.1007/s41365-022-00989-0 A novel 4D resolution imaging method for low and medium atomic number objects at the centimeter scale by coincidence detection technique of cosmic-ray muon and its secondary particles D.I.M. Setiadi , PSNR vs SSIM: imperceptibility quality assessment for image steganography . Multimed Tools Appl. 80 , 8423 - 8444 ( 2021 ). doi: 10.1007/s11042-020-10035-z http://doi.org/10.1007/s11042-020-10035-z PSNR vs SSIM: imperceptibility quality assessment for image steganography S. Agostinelli , J. Allison , K. Amako et al. , GEANT4 collaboration, GEANT4 — a simulation toolkit . Nucl. Instrum. Meth. A 506 , 250 - 303 ( 2003 ). doi: 10.1016/S0168-9002(03)01368-8 http://doi.org/10.1016/S0168-9002(03)01368-8 GEANT4 collaboration, GEANT4 — a simulation toolkit J. Feng , Z. Zhang , J. Liu et al. , A thermal bonding method for manufacturing Micromegas detectors . Nucl. Instru. and Meth. A 989 , 164958 ( 2020 ). doi: 10.1016/j.nima.2020.164958 http://doi.org/10.1016/j.nima.2020.164958 A thermal bonding method for manufacturing Micromegas detectors C. Hagmann , D. Lange , D. Wright , Cosmic-ray shower generator (CRY) for Monte Carlo transport codes . IEEE Nucl. Sci. Symp. Conf. Rec. 2 , 1143 - 1146 ( 2007 ). doi: 10.1109/NSSMIC.2007.4437209 http://doi.org/10.1109/NSSMIC.2007.4437209 Cosmic-ray shower generator (CRY) for Monte Carlo transport codes A.G. Bogdanov , H. Burkhardt , V.N. Ivanchenko et al . Geant4 simulation of production and interaction of Muons . IEEE Trans. Nucl. Sci. 53 , 513 - 519 ( 2006 ). doi: 10.1109/TNS.2006.872633 http://doi.org/10.1109/TNS.2006.872633 Geant4 simulation of production and interaction of Muons All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Si-Yuan Luo, Yu-He Huang, Xuan-Tao Ji, Lie He and Xiao-Dong Wang. The first draft of the manuscript was written by Si-Yuan Luo and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01070-6 This work was supported by the National Natural Science Foundation of China (No. 11875163) , Natural Science Foundation of Hunan Province (Nos. 2021JJ20006 and 2021JJ40444) , Ministry of Science and Technology of China (No. 2020YFE0202001) , and Department of Education of Hunan Province (Nos. 19B488 and 21A0281) . 2022-04-18 2022-05-25 2022-06-01 2022-07-12 2022-08-27 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01070-6 http://dx.doi.org/10.1007/s41365-022-01070-6 . 2022 Nuclear Science and Techniques 07 33