Basic principle of parallel computing

The primary concept of parallel computing is to divide a large task into subtasks for some processors and speed up the calculation; each processor works in concert with the others to conduct subtasks in parallel. MPI and OpenMP are primarily used in parallel computing for HPC to facilitate internode information delivery, resource allocation, and memory sharing.

Scalability, which refers to a system's ability to expand in response to future changes in requirements, is an important metric for evaluating parallel computing methods in HPC. Specifically, in this study, satisfactory scalability means that the parallel computing method is less time-consuming and ensures favorable efficiency when more projection images with more emitting particles are required. This study assessed the scalability and performance of the parallel computing mode by considering parallel speed-up and efficiency [21,22]. Speedup is defined as the ratio of the single-thread execution time to the execution time when one or more nodes with multiple threads are used. Parallel speedup, denoted as Sn, can be expressed as $\begin{array}{c}{S}_n=\frac{T_1}{T_n},\end{array}$ (2) where T₁ and T_n are the run time using one and n threads, respectively. Parallel efficiency, denoted as E_n, is defined as the average utilization of the n allocated threads and it can be expressed as $\begin{array}{c}{E}_n=\frac{S_n}{n},\end{array}$ (3) where n is the number of threads and Sn refers to the speedup with n threads.

Implementation of parallel computing

Parallel computing techniques have been widely used since the release of Geant4 (version 10.0). G4-MPI cross-node scaling and G4-MT cross-thread scaling are used to realize hybrid parallel computing supported by data collection methods such as tuples, scorers, and histograms. A master node and several slave nodes are used to implement the computing process. Each slave node is connected to several threads that cooperate to output the preliminary computing results, and once the results from each node are merged, they form the final output, which is sent to the master node. In this study, the final output is the energy deposition collected from the histograms.

In this study, a parallel interface and DR application were developed. The DR application in this study included a DR system and function-calling scripts that deployed parallel interfaces. In addition to incorporating the X-ray tube, flat panel detector, and test model in the geometric modeling, Monte Carlo simulation parameters such as particle definition, physical process, and data collection were specified by the imaging system to accurately simulate X-ray interactions with the imaged object. The parallel interface connected the parallel computing functions to a DR system by transmitting data to different nodes and threads. Relevant scripts were developed to construct the interface. A combination of nodes and threads was defined for the submission of the parallel input script, then the energy deposition was collected in the automated merging of the output data script.

Figure 2 shows the application of Geant4-based parallel computing to the simulation of a DR system. The blue boxes represent the DR system and parallel interface, whereas the yellow box displays the 3-D reconstruction process, all of which were developed by the methods described in this article. The red box represents the Geant4 kernel and the purple box represents the HPC developed by Chengdu HPC. The script extracts of the parallel input submission, energy deposition results collection, grayscale transform, and 3-D reconstruction are shown in Fig. 2 below.

Fig. 2Full-screen View

Overview of Geant4-based parallel computing applied in DR. Red pointer line: the DR system's parameters and scripts deploying the parallel interface are settled before the parallel computing. Pink pointer lines: the mapping relationship of parallel computing modes between the Geant4 kernel and HPC. Green pointer lines: the overall system operation flow which starts from the parallel interface deploying scripts and ends at the reconstructed 3-D image generation. The numbers in green indicate the sequences of main calculation and processing.

By employing the interface and application, the simulation of the DR was enabled by G4-MPI and G4-MT, which form the basis of the following sections.

Accuracy of parallel computing

This section verifies the accuracy of parallel computing by comparing the projection images obtained using parallel and non-parallel computing models. In the parallel computing mode, ten nodes with full threads, which is one case of the hybrid G4-MPI and G4-MT modes, were used, whereas in the non-parallel computing mode, only one thread was used in the simulation. 10 billion particles with an energy of 180 keV were emitted to ensure the convergence of the Monte Carlo simulation. Peak Signal to Noise Ratio (PSNR) and Root Mean Square Error (RMSE) [23,24] were utilized to measure the correctness of parallel images with non-parallel images as the benchmark. The RMSE can be calculated as $\begin{array}{c}\text{RMSE}\left(I_{\text{p}},I_{\text{np}}\right)=\sqrt{\frac{1}{m\times n}{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m\left(\left(I_{\text{p}}\left(i,j\right)-I_{\text{np}}\left(i,j\right)\right)\right)}}}^2},\end{array}$ (4) where I_p and I_np are the images obtained in the parallel and non-parallel computing modes, respectively. The PNSR can be calculated as $\begin{array}{c}\text{PNSR}=10\times {\text{log}}_{10}\frac{{\left({\text{MAX}}_I\right)}^2}{{\left(\text{RMSE}\right)}^2},\end{array}$ (5) where MAXI is the largest pixel value one image has.

Two images from two modes in the same projection angle are shown in Fig. 3, and it can be seen that the two images basically remain the same. The RMSE and PSNR between them are 3.83 and 36.5, respectively, indicating that the two images have a small difference in grayscale. This demonstrates that the results of the parallel computing mode based on Geant4 remain consistent with those of the non-parallel mode; thus, the following sections of the article can be presented.

Fig. 3Full-screen View

Projection images of two computing modes. (a) Projection image of parallel mode. (b) Projection image of non-parallel mode.

Performance evaluation of G4-MPI and G4-MT

Six test groups were designed using different combinations of test nodes and threads for G4-MPI and G4- MT in Chengdu HPC. To ensure the accuracy and repeatability of the simulation results, each test was repeated 30 times to obtain the population mean value μi and population standard deviation value δi [25]. $\begin{array}{c}{\mu }_i=\frac{1}{n}{\displaystyle \sum _{j=1}^{30}{x}_i}\left(j\right),\end{array}$ (6) $\begin{array}{c}{\delta }_i=\sqrt{\frac{1}{n-1}{\displaystyle \sum _{j=1}^{30}\left(x_i\left(j\right)-{\mu }_i\right)}},\end{array}$ (7) where xi(j) represents the j^th sample value (speedup or parallel efficiency) of test group i and n is the number of repetitions of each test. The confidence interval is calculated as in Eq. (8) [26]. $\begin{array}{c}\text{interval}\left(c_1,c_2\right)={\mu }_i\pm Z_{\alpha /2}\times {\delta }_i/\sqrt{n},\end{array}$ (8) where α is the significance level and c₁ and c₂ are the upper and lower bounds of the interval, respectively. The critical value, Z_α∕2, is calculated according to the normal distribution table considering 99.7% as the confidence level for each test [27]. The detailed results are shown in Table 2 and Table 3. The line graphs of parallel speedup and efficiency for six tests are shown in Fig. 4.

Fig. 4Full-screen View

Parallel speedup and efficiency of G4-MPI and G4-MT. (a) Parallel acceleration (S_n); (b) Parallel efficiency (E_n)

Figure 4 shows that the speedup increases linearly with the number of nodes or threads, whereas the parallel efficiency tends to stabilize at a relatively high value as the number of nodes or threads increases. Furthermore, according to Tables 2 and 3, with a 99.7% confidence level, the true values of the speedup and parallel efficiency of the six tests were within the confidence level for both G4-MT and G4-MPI. This demonstrates that both G4-MT and G4-MPI have good repeatability and efficiency in simulating the DR process. The speedup and efficiency of G4-MPI are marginally higher than those of G4-MT because G4-MPI adopts a shared memory mode, which improves computing by having each thread share information during Monte Carlo calculations.

Scalability evaluation of hybrid G4-MPI and G4-MT for DR

In this section, ten tests for the DR simulation are designed to evaluate the scalability of the hybrid G4-MPI and G4-MT modes. The G4-MT in this section ran in the full-thread mode for each computational node with 32 threads. The nodes of the 10 tests ranged from 1 to 10. Ten billion X-ray particles were emitted during each test. To illustrate the scalability evaluation, the runtime of one node using one thread was utilized as the benchmark time for 10 tests. The results of the parallel computing of the DR simulation are listed in Table 4, and the parallel speedup and efficiency of the ten tests are plotted in Fig. 5.

Fig. 5Full-screen View

(Color online) Scalability evaluation of hybrid G4-MPI and G4-MT. (a) Runtime and parallel speedup (S_n); (b) Parallel efficiency (E_n).

Table 7 and Fig. 5 demonstrate that as the number of nodes increases, the speedup deviates from linear growth and the parallel efficiency decreases within a certain range. According to Amdahl's Law [28], which describes how speedup increases theoretically when multiple processors are used in parallel computing, different nodes have parallel overhead time such as communication and waiting; therefore, when the number of nodes increases, parallel speedup will deviate from ideality, which leads to a decrease in parallel efficiency. This is acceptable in this study because the speedup can reach up to 151.36 and the reduction in runtime is significant. For example, in the next section, 180 projection images were required for reconstruction. Based on the runtime required for one projection image using one thread, it can be estimated that 2,323.4 h are required. However, using G4-MPI in full threads, the process will only cost approximately 15.35 h in total. The strong scalability of G4-MPI in the full-thread mode for DR simulation was proven.

3-D reconstruction using filtered back projection algorithm

The FBP algorithm is widely used. It is a spatial processing technique based on Fourier transform theory [29]. One of its characteristics is that the projection under each acquisition projection angle is convolved by specific filters before back-projection, such that the toroidal artifact caused by the point spread function can be reduced. This study adopted the FBP algorithm to obtain the 3-D reconstructed images of the model. The FBP process is described as follows.

• Step 1: Obtain one 300 × 180 original projection matrix using the same rows of pixel information in all the original projection images.

• Step 2: Pre-process the original projection matrix data for calibration.

• Step 3: Conduct a one-dimensional Fourier transform for the result of Step 2.

• Step 4: Perform convolution filtering on the result of Step 3.

• Step 5: Conduct a one-dimensional inverse Fourier transform for the result of Step 4.

• Step 6: Perform a direct inverse projection on the result of Step 5 for the 2-D reconstructed slicing image.

• Step 7: Repeat Steps 1 to 6 300 times to produce 300 2-D slicing images.

• Step 8: Utilize the 300 images to obtain 3-D reconstructed images.

Step 1 considers the Nth row of the pixel information of each original projection image. N is initially equal to one and increases by one with each repetition. Thus, 300 original 300 × 180 projection matrices are generated after 300 repetitions. Three hundred 3-D reconstructed slicing images, each 300 × 300 pixels in size, are generated.

As the toroidal artifact is one of the main factors affecting the quality of 3-D reconstructed images, this study adopted a projection data pre-calibration method combining polynomial fitting and the probability statistics correction method in Step 2 [30]. The range of alternative correction factors was determined by polynomial fitting of the projected data column-by-column, and correction factors were determined using the maximum probability principle.

Among the eight steps above, the design of the filter in Step 4 has a significant impact on the final reconstruction results [31]. The selection of filters is actually the selection of window functions and to gain better-reconstructed image resolution, the inverse Fourier transform of the window functions should have a high and narrow central protrusion [32]. To obtain reconstructed images with better quality, this study used the three most commonly used filters: the Ram-Lak, Shepp-Logan, and Hamming filters to perform convolution filtering on the result of the one-dimensional Fourier transform. The 2-D reconstructed slicing image under 0° rotation is consistent with that shown in Fig. 1.

Quality evaluation of reconstructed images

To quantitatively measure the quality of the reconstructed images, three commonly used image quality assessment metrics, contrast-to-noise ratio (CNR) [33], average gradient (AG) [34], and image entropy (IE) [35], were used. CNR is defined as the ratio of peak signal intensity to background intensity and is an objective indicator of average image quality. Its value is proportional to the image quality, and a higher CNR value indicates a better recognition of target defects. AG is sensitive to the ability of an image to express small details in contrast. In a certain direction of an image, AG tends to be larger with a greater change in gray level. Thus, AG could be used to measure the clarity of an image and reflects small detail contrasts and texture transformation features in an image. IE reflects the amount of average information in the image and represents the aggregation characteristics of the image grayscale distribution. The greater the entropy of the image, the richer the pixel grayscale contained in the image, and the more uniform the grayscale distribution.

The results of the reconstructed image quality assessment of the model used in this article are shown in Table 5. Based on the spatial characteristics of the model, this study presented the results of a quality assessment of the reconstructed images of the 100^th, 150^th, and 200^th layers of the test model. Reconstructed slicing images of the three layers depict the representation of hollow defects and solid aluminum. The quality of the reconstructed images varied considerably for different filters and the number of slicing layers in the reconstructed images.

First, the reconstructed slicing images using the Ram-Lak filter had better quality than those using the Shepp-Logan or Hamming filters. Second, there was a high degree of similarity between the metric values of the 100^th and 200^th layers because of the consistent physical structures of the two slices in the model, such as the density, material, and shape of the defect. Moreover, the CNRS of the 100^th and 200^th slices are smaller than that of the 150^th slice, mainly because the defect area in the 150^th slice of the image is smaller than those of 100^th and 150^th slices. This demonstrates that the 3-D reconstruction method based on parallel computing has good recognition of small defects. Furthermore, the AG and IE of the 100^th and 200^th slices were larger than those of the 150^th because of the greater change in the grayscale value of the former. It can be concluded that the images of the 100^th and 200^th slices showed more detail in contrast. In general, the value of the three metrics meets the desired expectation for the designated model. Figure 7 shows the reconstructed 2-D images using the three aforementioned filters. It is evident that the toroidal artifacts are more prominent in the images reconstructed with the Shepp-Logan and Hamming filters compared to those reconstructed with the Ram-Lak filter.

Fig. 7Full-screen View

Representative reconstructed slicing images a1-a3: Reconstructed images of 100^th slice b1-b3: Reconstructed images of 150^th slice c1-c3: Reconstructed images of 200^th slice a1,b1,c1: Reconstructed images of Hamming filter a2,b2,c2: Reconstructed images of Shepp-Logan filter a3,b3,c3: Reconstructed images of Ram-Lak filter

$\text{CNR}=\frac{\left|{\mu }_i-{\mu }_o\right|}{\sqrt{{\delta }_i{}^2+{\delta }_o{}^2}}$, ${\delta }_i^2=E{\left\{\left(|s_i{|}^2-{\mu }_i\right)\right\}}^2$, ${\delta }_o^2=E{\left\{\left(|s_o{|}^2-{\mu }_o\right)\right\}}^2$, μi and μo are the average pixel value of the region of the interest (ROI) and background, respectively, and ${\delta }_i$ and ${\delta }_o$ are the variances of ROI and background, respectively. $\text{AG}=\frac{1}{M\times N}{\displaystyle \sum _{i=1}^M{\displaystyle \sum _{j=1}^N\sqrt{{\left(\frac{\partial f\left(i,j\right)}{\partial x}\right)}^2+{\left(\frac{\partial f\left(i,j\right)}{\partial y}\right)}^2}}}$, whereas M×N is the pixel size of the image. $\frac{\partial f\left(i,j\right)}{\partial x}$ and $\frac{\partial f\left(i,j\right)}{\partial y}$ are the differences of the pixel value in the x and y directions, respectively. $\text{IE}=-{\displaystyle {\sum }_{g=0}^{L-1}{P}_g\times {\text{log}}_2{P}_g}$, L is the overall grayscale level of the image, and $P_g$ is the probability of gray level g.

To further discuss the quality of the 3-D-reconstructed images, we compared the defect situation of the model in the ideal case with that of the reconstructed slicing image sets. As the test model was accurately modeled using Geant4, this study utilized the parameter table as an inspection standard for the set of 300 reconstructed slicing images.

Among the 300 images were some blank images because the size of the flat-panel detector was larger than that of the test model. During this process, the X-ray particles were captured by a flat-panel detector without decay. In this case, after the data of the associated units on the flat-panel detector were processed, the corresponding pixels appeared blank.

In this process, excluding 100 blank images, the model was divided into 200 slices. A total of 148 slices contained small square and cylindrical defects, whereas 52 slices contained only cylindrical defects. The ratio of layers containing the square defect and non-blank layers (r₁) is 0.74:1, whereas the corresponding ratio (R₁) according to parameter Table 1 is 0.75:1. The ratio of the length of the side of the small square defect to the length of the side of the large square (r₂) in Fig. 7 is 0.37489:1, whereas the corresponding ratio (R₂) according to parameter Table 1 is 0.375:1. The ratio of the diameter of the round defect to the length of the side of the large square is 0.0749:1 (r₃), whereas the corresponding ratio according to parameter Table 1 is 0.075:1 (R₃). The difference between the three ratios r₁, r₂, and r₃ is within the acceptable error range. Examining the images of the remaining layers showed that the results were similar (Table 6).

Fig. 8 illustrates the grayscale value variation curves of the same row of pixels of the two layers of the reconstructed slicing images. These two images were selected according to the parameter table whose defect situation was theoretically consistent, such as two adjacent layers. Owing to the hollow defects inside the model, the reconstructed sliced images exhibited obvious differences in the grayscale values for specific row pixels. From Fig. 8, it can be concluded that the grayscale value varied between [0,1] with the appearance of defects in the image, and that of the reconstructed slicing image coincided closely with the adjacent layer of the reconstructed slicing image. The RMSE values of the three sets of pixel data were 0.000378, 0.000230, and 0.000496, indicating an outstanding match within each set of pixel data.

Fig. 8Full-screen View

Grayscale value variation curves. (a) 200^th row pixel position; (b) 150^th row pixel position; (c) 100^th row pixel position

To obtain the complete 3-D reconstruction view shown in Fig. 9, a 3-D model was formed by stacking 180 2-D slice images in the correct order in 3-D space. However, inconsistencies in the images can result in a blurred appearance of the final 3-D model. To address this issue, a smooth filtering operation was applied to obtain a highly restored original model. The reconstructed surface geometry had almost no holes, indicating a high level of accuracy in the 3-D reconstruction process. It can be seen that the reconstructed model reproduced the defect features and external shapes modeled in Geant4 very well. The disadvantage of this method is that some irregularities exist on the surface, which may be caused by incomplete elimination of toroidal artifacts.

Fig. 9Full-screen View

Reconstructed image of 3-D view of the test model