Hi'CT: a pixel sensor-based device for ion tomography

Carbon ions, commonly referred to as particle therapy, have become increasingly popular in the last decade. Accurately predicting the range of ions in tissues is important for the precise delivery of doses in heavy-ion radiotherapy. Range uncertainty is currently the largest contributor to dose uncertainty in normal tissues, leading to the use of safety margins in treatment planning. One potential method is the direct relative stopping measurement (RSP) with ions. Heavy-ion CT (Hi′CT), a compact segmented full digital tomography detector using monolithic active pixel sensors, was designed and evaluated using a 430 MeV/u high-energy carbon ion pencil beam in Geant4. The precise position of the individual carbon ion track can be recorded and reconstructed using a 30 μm × 30 μm small pixel pitch size. Two types of customized image reconstruction algorithms were developed, and their performances were evaluated using three different modules of CATPHAN 600-series phantoms. The RSP measurement accuracy of the tracking algorithm for different types of materials in the CTP404 module was less than 1%. In terms of spatial resolution, the tracking algorithm could achieve a 20% modulation transfer function normalization value of CTP528 imaging results at 5 lp/cm, which is significantly better than that of the fast imaging algorithm (3 lp/cm). The density resolution obtained using the tracking algorithm of the customized CTP515 was approximately 10.5%. In conclusion, a compact digital Hi'CT system was designed, and its nominal performance was evaluated in a simulation. The RSP resolution and image quality provide potential feasibility for scanning most parts of an adult body or pediatric patient, particularly for head and neck tumor treatment.


Introduction
Heavy-ion and proton therapies are considered the most advanced particle therapy technologies. As research and clinical applications progress, they have become reliable, safe, and promising treatment options. Wilhelm Conrad Röntgen first discovered X-rays in 1895 [1], and Robert R. Wilson proposed accelerating charged particles, specifically protons, for cancer treatment in 1946 [2]. In 1954, the first patient was treated with proton therapy at the University of California, Berkley [3]. In 1990, the first hospital-based proton therapy center was opened at Loma Linda University [4]. At the Lawrence Berkeley Laboratory (LBL), Cornelius A. Tobias proposed that the use of heavy ions improves relative biological effectiveness (RBE). In 1975-1992, LBL explored the possibility of using various types of charged particles, including He, N, O, C, Ne, Si, and Ar, for treatment, and carbon ions were discovered to be the best option [5][6][7]. By the end of 2021, more than 300,000 patients have received particle therapy, with the equipment for the said therapy also improving rapidly [8][9][10].
China has been investigating heavy-ion therapy for the past two decades. Since then, the development of particle therapy devices has accelerated in China. The Institute of Modern Physics (IMP) of the Chinese Academy of Sciences, which plays a leadership role in this area, established the first Heavy Ion Research Facility in Lanzhou (HIRFL) in 1989. The carbon-ion beams provided by HIRFL were used 111 Page 2 of 11 in cell and animal experiments in 1993, and the possibility of clinical application was explored in 1995. In 2006, a superficial tumor terminal for heavy-ion therapy was successfully developed in collaboration with hospitals for the clinical application of human superficial malignancies. In 2008, the IMP completed the cooling storage synchrotron (CSR) at the Lanzhou Heavy Ion Accelerator, and a deep tumor treatment terminal was subsequently developed and used in clinical trials. Between 2006 and 2013, more than 200 patients underwent carbon-ion radiotherapy (CIRT) at the HIRFL. Since then, China has become the fourth country, after Japan, Germany, and the USA, to initiate research on heavy-ion therapy. Based on the results, IMP built the heavy-ion medical machine (HIMM) in Wuwei and Lanzhou City, Gansu Province in 2012. The device was successfully launched in Wuwei in 2015, and clinical trials with 47 patients were completed in 2019. It was officially adopted for clinical application in March 2020. Three-hundred patients have been treated in HIMM, and good therapeutic effects and tolerability were observed until 2021. Since then, particle therapy devices in China have developed rapidly [11,12].
When creating a treatment plan, the range of heavy ions in the tissues plays an important role in successfully delivering an accurate dose to the target volume. During this process, the X-ray CT images, which must be converted into relative stopping power (RSP) of heavy ions, are normally used in heavy-ion treatment planning. The conversion process leads to uncertainties of approximately 2% in the further range prediction of heavy ions, which corresponds to an absolute range uncertainty of 2 mm at a depth of 10 cm in the target volume [13,14].
These range uncertainties are clinically handled by irradiating a "safety margin," which is an extended margin around the tumor, to ensure full coverage of the target volume. The margin introduces a higher probability of excessive doses to normal tissues. This risk can be reduced by avoiding beam control points where normal tissues can be irradiated; nevertheless, fewer degrees of freedom of the beams can result in a potential full prescription dose [15], which can lead to severe radiation side effects.
Several techniques have been suggested to reduce range uncertainties [16]. They can generally be classified into one of two categories: in vivo dose measurement, such as prompt gamma emissions, and direct stopping power measurement, which is the focus of this study. Instead of converting X-ray CT to heavy-ion RSP, direct measurement of heavy-ion range or energy loss provides a potential feasibility of precisely predicting the heavy-ion range in the human body. Heavier ions have fewer multiple scattering effects owing to their larger momentum, making the straight-line approximation easier to implement. This allows the use of standard tomography reconstruction methods from X-ray CT.
With the progress of modern detector technologies that are mostly developed for particle physics, several types of compact detection systems based on scintillation fibers [17], GEM [18], SiPMs [19], and CMOS APS sensors [20,21] have been investigated by tracking each particle and measuring its remaining energy or range [22,23].
Ion CT requires an accelerator and a gantry to deliver the beam for tomography imaging. Because these are extremely expensive, even in the context of medical scanners, the only feasible approach is to develop a compact system that can be used in existing clinical particle therapy facilities. Clinical heavy-ion accelerators for therapy typically exhibit a maximum energy of approximately 430 MeV/u, corresponding to a range of approximately 31 cm in water. This would be sufficient for scanning most parts of an adult body or pediatric patient, particularly the head and neck; however, it cannot be scanned from all angles at this energy.
In this study, a compact digital heavy-ion CT (Hi'CT) system was designed, and its performance was evaluated using a Geant4-based simulation. The digital detector consists of a tracker and a range-measurement sampling calorimeter, which are fully based on monolithic active pixel sensors. Using a stack of range detectors with micron-level positionresolution pixel sensors, accurate carbon ion tracks can be recorded and reconstructed by the system. The detector performance and several standard phantoms contained in CAT-PHAN 600 series, which are typically used for X-ray CT, were tested to verify its RSP accuracy and imaging qualities.

Detector conceptual design
The detector modeled in the simulation consists of two key components: a multilayer tracker used for reconstructing the angle and a sampling calorimeter for recording the residual range of carbon ions after penetrating the imaging object. Four tracking layers, which are divided into two groups, are located before and after the imaging object. A 10-mm-thick tungsten block, which has a large stopping power, is placed between the tracker and calorimeter to decrease the energy of ions. Twenty sampling layers are placed in the calorimeter immediately after the tracker, each of which is composed of absorbers, silicon sensors, printed circuit boards (PCB), and air gaps. Monolithic active pixel sensor (MAPS), a unique technology for high-granularity detectors in particle physics, is used for the sensitive component of the detector, which contains 640 × 640 pixels with a squared pixel size of 30 μ m. In each layer, 40 MAPS sensors are mounted next to each other in two rows, resulting in a 0.01-mm narrow gap between the sensors along the x and y directions (Fig. 1).
Detector performance and budget are the main considerations in the conceptual design. The sampling calorimeter aims to cover the clinical ion energy that passes through a certain thickness (e.g., 10-25 cm water-equivalent thickness (WET) for the head and neck) of the human body. Therefore, optimization of the number of sampling layers and the thickness of the absorbers in the calorimeter is required. The parameters of the entire system should match the following expressions: where In the above expressions, T in is the total WET of the incident carbon beam, and X up and X down are the upper and lower WET limits of the object to be measured, respectively. T ab represents the WET of the single-layer absorber, whereas T det is the WET of the single-layer detector. W, D and n are the WET of the tungsten block, the WET of sampling layers, and the number of sampling layers, respectively.
Parameters for clinical applications must meet specific physical conditions. The Hi'CT system is based on the logic that ions can penetrate the patient's body, and the Bragg peak can eventually be stopped in the calorimeter layers after the tungsten block. The number of layers n must satisfy the condition of Eq. 3.
W and T ab should satisfy Eqns. 4 and 5. (1) When T in is known from the incident carbon ions, the X up and X down are calculated according to the typical WET of the object to be imaged. T det can be obtained based on the pixel sensor parameters, for example, the thickness of the epitaxial layer and materials. After finding reasonable thickness tungsten block (W ) via simulation, the possible optimization values for n and T ab of the detector can be easily determined. These theoretical equations were used to reduce the time required for exhaustive exploration. It is also helpful in estimating the cost of the final prototype design.
In the current design, factors such as detection accuracy, cost, and production process were comprehensively considered. The system parameters obtained based on the selection formula are listed in Table 1.

Monte-Carlo simulation
The detector was implemented using the Monte-Carlo software toolkit Geant4 [24] with the FTFP_BERT physics list. The simulation and reconstructed analysis software were customized and run on a computer platform with an Intel Xeon Gold 6248R CPU and 128 GB RAM. The energy deposit in each sensor plane was discretized into a list of fired pixels analogous to real fire pixels. An energy threshold of 10 MeV was used to simulate sensor sensitivity controlled by the reference voltages.
The beam characteristics in the simulation were based on the parameters of the Wuwei Carbon Ion Therapy System. An initial 430-MeV/u carbon pencil beam was simulated. The position spread assumed a Gaussian distribution with a full width at half maximum (FWHM) of 3 mm. Figure 2 shows the geometric information of the detector setup in Geant4. The beam, which is shot perpendicular to the upstream detector, was in the scan mode along the horizontal plane with a step length of 0.1 mm. Owing to limited data processing power, 100 carbon ions were generated in every step, which corresponded to a central fluence of

Water equivalent path length calibration
To verify the direct relationship between RSP from the Hi'CT system and the water-equivalent path length (WEPL), the WET values of the detector were calibrated in the simulation. The calculations were performed using the same carbon-ion beam parameters described in Sect. 2.2. By calculating the ratio in Eq. 6 [25] between the thickness of the water phantom and the equivalent stop depth of a given material in which 1000 carbon ions are fully stopped, the reference RSP of a given material can eventually be defined.
where T w and T m are the thicknesses of the water phantom and a given material, respectively.

Tracking and image reconstruction algorithm
In this study, several algorithms were developed to transfer pixel response hits into the carbon tracks used for the final image reconstruction. The main goal of the tracking algorithm was to provide sufficient precision for the carbon tracks before and after penetrating the phantom with the recorded matter. The energy deposition of Geant4 steps within a pixel was summed up, and this is called a "hit". A 10-MeV energy cut of secondary particles was applied to simulate the sensitivity of pixel sensors, which is controlled by the reference voltage in practice. Two different reconstruction algorithms-tracking and nontracking-were developed to verify the reconstructed image qualities for further clinical applications. The results were compared with pure carbon tracks selected in Geant4, which was treated as an ideal case.
The tracking algorithm provides more precise carbon track accuracy. However, the time consumption and relatively low efficiency are the main limitations of the current stage. Therefore, an algorithm called the "fast" imaging algorithm is developed for prescanning the imaging object without tracking the intermediate step applied.
From these events, selected tracks were reconstructed according to the following procedure: • For every single carbon event, all the hits it produced in the detector with energy above 10 MeV are selected as candidates. • The initial "seed" position of the carbon event is determined in this step according to the (T, ) of the beam, where T and are the offset and rotation angle of the phantom, respectively. • Use the pixel whose Euclidean distance r is the least as the seed for the next layer. • For a given layer, assign all the hits within r = 1 mm into the tracklist until no hit is found in the next two consecutive layers.  Several conditions must be clarified during the tracking process. If the beam does not pass through the tungsten block, the endpoint is defined as layer 0. If there is only one layer without a candidate pixel hit between consecutive layers with candidates, the tracking process continues until reaching the next layer. In some cases, ions reach the maximum depth of the detector, which frequently occurs around the edge of the imaging object; the range of the track is marked by the maximum number of layers of the detector. Statistical effects must also be considered when obtaining the WEPL values at each coordinate. By plotting the histogram distribution, the range of all particles at a single coordinate can be determined. There will be a peak position near the end of the beam. It provides an excellent idea of where the Bragg peak position is (i.e., the number of layers in the corresponding sampling layer). By setting a threshold for the histogram distribution, information about the peak position can be obtained, which can be either the response of one or two sampling layers. If there is only one response, the WEPL is calculated in the middle of the silicon sensor corresponding to the sampled layer. If there are two responses, the average of the WEPL calculated at the middle of the silicon sensors corresponding to the two sampled layers is chosen as the WEPL at that coordinate.
In contrast to the tracking algorithm, the fast mode aims to reconstruct the image by directly using information from events, which significantly increases data processing efficiency and can be used for fast diagnostic purposes. It reconstructs the sinogram according to the following procedure: • For every single carbon event, all the hits produced in the detector with energy above 10 MeV are selected as candidates. • The initial "seed" position of the carbon event is determined according to the (T, ) of the beam. • The histogram distribution of the number of hits in each layer at each coordinate is obtained. • The layer where the peak point of the histogram is located is used as the endpoint of the carbon beam range at this coordinate.
The preprocessing of the data between the Monte-Carlo simulation and the reconstruction was performed using the Python Pandas library. Parallel computing enables the reconstruction process to be executed within a reasonable amount of time. In this study, the mapping from detector space to sinogram space was performed in a way that yields the best image quality. Because the sinogram space is a grayscale image, each coordinate in the sinogram space must be filled with 1-dimensional data using the beam range information, that is, WEPL. To transform the detector space into sinogram space, the data at the event and coordinate levels must be processed separately. All images under these three scenarios were reconstructed using the filtered back projection (FBP) algorithm with the Ramachandran-Lakshminarayan (Ram-Lak) filter.

Phantom imaging
The performance of the simulated Hi'CT scanner was validated using CTP404 and CTP528 from the CATPHAN 600 phantom (The Phantom Laboratory, Salem, NY, USA). All these phantoms were modeled in SolidWorks and transferred with CADMesh version 2.0.3, a CAD interface, to Geant4 [26]. The CTP404 module ( Φ = 15 cm, Height = 2.5 cm) consists of eight cylindrical cavities with a diameter of 1.3 cm. The cavities are filled with different materials such as Teflon, Delrin, acrylic, polystyrene (PS), polymethylpentene (PMP), low-density polyethylene (LDPE), air, and water ( Fig. 3 (top left)). The base material of the CTP404 module is epoxy. The RSP evaluation of the reconstructed images was conducted using ImageJ version 1.46r, a Javabased open-source image-analysis software [27]. The insert materials of the CTP515 module ( Φ = 15 cm, Height = 4 cm) are muscle, eye lens, and bone, and the base material is adipose. The insert materials of the CTP515 module were modified to explore head and neck imaging possibilities. The insert material of the CTP528 module ( Φ = 15 cm, Height = 4 cm) is aluminum, and the base material is epoxy. The CTP528 line pair module provides 1 through 21 aluminum line pairs per centimeter, such that all patterns share the same distance from the center of the phantom. This module was used to evaluate the spatial resolution measured in the phantom images reconstructed from the simulated data. For a quantitative comparison of the spatial resolutions, the modulation transfer function (MTF) was calculated for both cases using a custom script written in Python. For each of the line pair patterns with lp = 1-5, the average reconstructed RSP of the aluminum peaks and the average RSP of the base material troughs were calculated. The MTF for each analyzed line-pair group was calculated using the following equation [28]: where ⟨ r p ⟩ and ⟨r t ⟩ are the mean RSP values of the aluminum and base within the region of a given line pair, respectively, and ⟨r Al ⟩ and ⟨r b ⟩ are the RSP values of the aluminum and epoxy base measured by the water phantom in the simulation, respectively. It is worth mentioning that for the given region of a line pair, a rectangular area of approximately 80% of the size of a line pair was selected within the boundaries of each insert, and the mean RSP and standard deviation were calculated using the standard ImageJ tool.
The CTP515 module consists of a series of cylindrical rods of various diameters and three contrast levels to measure low-contrast performance. The 40-mm-long rods provide consistent contrast values at all z-axis positions, thereby avoiding volume-averaging errors. Because the original contrast of the rods (0.3, 0.5, and 1.0%) was significantly small to be visualized in the Hi'CT system, the materials in the rods were replaced by three different types of materials similar to the human body, resulting in density contrasts of 10.5, 12.6, and 94.7%, respectively.

Tracking accuracy
Imaging reconstruction speed and quality are sensitive to the accuracy and stability of the tracking process. The imaging dataset of CTP404 at 0 • was chosen as the candidate dataset to verify the tracking accuracy. Figure 4 shows the difference in the number of stopping layers between the tracking and ideal cases. The good agreement of the calculation resulting in the accuracy of the tracking is more than 94%. Figure 5 shows the sinogram and the reconstructed image of the CTP404 module, which were used to evaluate the profile of the module. As shown in the image, eight cylinders on the outer lane side and four rods are visible for the three different scenarios, whereas acrylic spheres with various diameters close to the center are indistinguishable. For image quality, particularly the edge sharpness of the cylinders, the image with the tracking algorithm is slightly  Fig. 4 Histogram of the range difference between the ideal and tracking cases less than the ideal pure carbon track reconstruction but better than that of the fast-imaging algorithm. In all three cases, evident annular and strip-shaped artifacts can be observed, which are mainly because of using the center of the layer as the range position during the range calculation. These artifacts can also be observed in other modules (e.g., CTP528 and CTP515).

RSP resolution
Although spatial resolution may have the greatest impact on how the image appearance, RSP resolution and accuracy in an Hi'CT system are more important for predicting the ion range for treatment planning. To calculate the RSP values of various materials, the diameter of the region of interest (ROI) in the reconstructed image was constrained to 64 pixels (corresponding to 11.38 mm on the image scale), which accounted for approximately 8% of the phantom diameter (150 mm). Table 3 lists the values of RSP measurements of the cylinders and rods in the CATPHAN 404 module. Figure 5 shows three reconstructed images of the CTP404 module using the ideal, tracking, and fast imaging algorithms (no-tracking). From the image, it can be concluded that eight peripheral and four internal pillars can be clearly reconstructed; however, the internal low-density globules cannot be imaged. From the perspective of the structural integrity of the interpolator, the image quality of the tracking algorithm is significantly improved compared with the fast-imaging method and is close to the ideal situation. However, regarding the uniformity of the interpolated image, the fast algorithm is worse than the tracking algorithm on some high-density interpolated objects, such as Teflon.
The RSP resolution of the eight interpolations was evaluated by defining a circular ROI that covers a diameter of approximately 8% of the image size. The results are shown in Table 3. The RSP resolution of the tracking algorithm can be lower than 1%. Although the effect of  the fast algorithm is slightly worse than that of the tracking algorithm, it benefits from the improvement of range measurement accuracy owing to the thin-layer design of the silicon pixel detector, which makes the RSP resolution of interpolators less than 2%, except for PMP. Figure 5 shows the reconstructed images of the CTP404 module using the simulated data. Table 3 lists the reconstructed RSP values for the insertions into the module. The relative difference between the ideal and tracking RSP is below 1% for all materials, whereas the Fast RSP is below 2% for all materials except PMP. The difference between RSP ref and Fast RSP for PMP is 2.43%, and the reconstructed RSP is in relatively good agreement with the expected value.

Spatial resolution
The spatial resolution of the system was evaluated using the CTP528 module. Figure 6 shows the carbon-ion tomography image results of the CTP528 module under three conditions: ideal, tracking, and fast imaging algorithms. It can be observed that the line-pair structure in the module can be completely reconstructed in all three cases. Figure 7 shows the MTF curve calculated for the first five groups of line pairs; the results of which agree well with the ideal situation, and both reach 20% at 5 lp/cm. The fast-imaging algorithm reaches approximately 10.5% points only at 3 lp/cm, owing to the range uncertainties caused by the interference of secondary particles.
The customized CTP515 phantom was used to evaluate density resolution. To simulate the main tissue materials of the human head and neck, the base material of CTP515 is made of adipose tissue, and the three materials with different densities are muscle, lens, and bone. Figure 8 shows the carbon-ion tomography image results of the CTP515 module under three conditions: ideal, tracking algorithm, and fast imaging. From a visual point of view, these three methods can better distinguish high-contrast tissues, such as bone and soft tissue. As the number of cases decreases, the tracking method exhibits more serious artifacts than the ideal case. The tracking algorithm can distinguish between muscles, adipose, and lens with similar densities to a certain extent. However, the fast algorithm cannot effectively reconstruct the geometric structure and image contrast.

Discussion
According to the results, both the tracking and fast algorithms can perform complete imaging of the internal structures of the CTP404 and CTP528 modules. As an accurate algorithm, the tracking algorithm yields better imaging results than the fast algorithm. For RSP resolution, the RSP measurement accuracy of the tracking algorithm for different types of materials in the CTP404 phantom is less than 1%, which meets the clinical requirements for distance-to-agreement (DTA)<1% in proton heavy-ion imaging [29]. In terms of spatial resolution, because the secondary particles in the fast-imaging algorithm increase the uncertainty of the algorithm, the tracking algorithm was used to achieve a 20% normalization MTF value of CTP528 imaging results at 5 lp/cm, which is significantly better than the spatial resolution of the fast-imaging algorithm (10.5% at 3 lp/cm). To verify the imaging effect of the detection system on the main materials of human tissues, the CTP515 custom phantom composed of the main tissue components of the head and neck was used to evaluate the density resolution, which is approximately 10.5% for the tracking algorithm. On the one hand, owing to the influence of artifacts and noise, the edge sharpness of the inserts in the module decreases. On the other hand, owing to the low number of effective cases after case selection, the detection ability of low contrast is reduced. However, it is worth mentioning that this study used only the most basic FBP algorithm with GPU acceleration to evaluate the intrinsic performance of the images. The impact of different types of algorithms on the image reconstruction quality will be discussed in future work.
Fast mode bypassed the tracking process and used the peak position of the raw-fired pixel hit count distribution to approximately locate the peak position of the Bragg peak. As a result, This algorithm led to a significant reduction in the reconstruction speed by a factor of 10, whereas the error remained within acceptable limits. The fast mode is important for clinical applications and can be used to prescan patients. Because the fast mode is based on a different idea than the tracking algorithm, the optimization of the algorithm is also different from that of the tracking algorithm, which will be discussed in future work.
From the reconstruction perspective, the current reconstruction process may still be two orders of magnitude slower than X-ray reconstructions because of the low efficiency of the data preselection and tracking process, which is still far from the clinical requirement of 5 min. However, this problem can be solved by optimizing the tracking and reconstruction algorithms, and increasing the computing power in parallel. Moreover, introducing an iterative reconstruction algorithm can further reduce the required number of scanning angles, reduce artifacts, and potentially decrease the dose to the patient. Overall, the results presented in this study require experimental validation, and further improvements are likely to be made based on these experimental results.
The quality of the reconstructed image is affected by strip and ring artifacts. Annular artifacts are mainly caused by multiple strip artifacts, whereas strip artifacts may be mainly because of the paragraph type detection structure causing the discreting of range definition. Although artifacts do not have a significant impact on the differentiation of the main structures in the phantom, they lead to a decrease in the RSP evenness and edge sharpness for different structures in the phantom. Reducing the influence of artifacts on imaging quality will be discussed in future research.

Conclusion
In this study, a compact digital heavy-ion tomography system Hi'CT based on multilayer MAPS pixel sensors was proposed and evaluated under a 430-MeV/u clinical carbon-ion pencil beam in Geant4. Consistent with past literature, this study demonstrated that the Hi'CT system can be designed with high stopping power accuracy. Benefitting from the 30 μ m × 30 μ m small pixel pitch size, the detector can accurately reconstruct the carbon-ion tracks before and after penetrating the object to be imaged, which is used to solve the range uncertainty problem caused by the use of traditional X-ray CT in heavy-ion planning. In the simulation, carbonion tracking and fast algorithm were established. Two standard modules (CTP404 and CTP528) from the CATPHAN 600 series and a customized CTP515 module using human tissue-like materials were used for tomography based on the FBP method and compared with the ideal situation of eliminating secondary particle interference.