Beam Angle Optimization for Double-Scattering Proton Delivery Technique Using an Eclipse Application Programming Interface and Convolutional Neural Network

To automatically identify optimal beam angles for proton therapy configured with the double-scattering delivery technique, a beam angle optimization method based on a convolutional neural network (BAODS-Net) is proposed. Fifty liver plans were used for training in BAODS-Net. To generate a sequence of input data, 25 rays on the eye view of the beam were determined per angle. Each ray collects nine features, including the normalized Hounsfield unit and the position information of eight structures per 2° of gantry angle. The outputs are a set of beam angle ranking scores (S beam) ranging from 0° to 359°, with a step size of 1°. Based on these input and output designs, BAODS-Net consists of eight convolution layers and four fully connected layers. To evaluate the plan qualities of deep-learning, equi-spaced, and clinical plans, we compared the performances of three types of loss functions and performed K-fold cross-validation (K = 5). For statistical analysis, the volumes V27Gy and V30Gy as well as the mean, minimum, and maximum doses were calculated for organs-at-risk by using a paired-samples t-test. As a result, smooth-L1 loss showed the best optimization performance. At the end of the training procedure, the mean squared errors between the reference and predicted S beam were 0.031, 0.011, and 0.004 for L1, L2, and smooth-L1 loss, respectively. In terms of the plan quality, statistically, PlanBAO has no significant difference from PlanClinic (P >.05). In our test, a deep-learning based beam angle optimization method for proton double-scattering treatments was developed and verified. Using Eclipse API and BAODS-Net, a plan with clinically acceptable quality was created within 5 min.

Recently, with increasing interest in proton and heavy ion therapy, which rely on the characteristics of a Bragg peak and a relatively higher radiation biological effect than X-ray therapy, several recent studies on intensity-modulated ion therapy (35) in BAO research have been published (1,(36)(37)(38)(39).
The present study was inspired by two previous studies on BAO. In 1999, Hosseini-Ashrafi et al. conducted a study on the BAO of X-ray therapy in which they used an artificial neural network (ANN) (30). In that study, the radiation treatment plans were divided into several templates, and the ANN classified the test data according to the template. The ANN consisted of three layers of a multi-layer perceptron. The input contained 12 precalculated features, which were the body contour outline, treatment volume, sensitive organs, and border of tissue inhomogeneity for each case. The output contained three types of binary data for eight classification tasks. The ANN was validated using the leave-one-out method, which showed the feasibility of applying ANNs to the BAO problem.
In 2002, Pugachev et al. published a research paper on BAO for IMRT (4). They proposed beam's eye view dosimetrics (BEVD) to overcome low computation speed, which is a disadvantage of the simulated algorithm. The BEVD score was calculated by using the geometric and dosimetric information of the patient. This score was used for ranking information and as a prescreening tool to optimize the beam orientation by using a simulated annealing-based BAO algorithm. The treatment plans generated with the guidance of the BEVD score were compared with those created with five equiangular-spaced beams. They validated the feasibility of the BEVD score for the BAO problem. The BEVD guidance indicated that the computational efficiency increased by a factor of~10.
In the current study, we developed a deep-learning based BAO method for the three-ports proton double-scattering (DS) technique using the geometric information of the patient computed tomography (CT) anatomy and Hounsfield unit (HU) data as well as a convolutional neural network (CNN).
In particular, the DS technique was used for large-field proton therapy. The proposed method requires only geometric information without a fluence optimization process. The geometric information is automatically extracted using an application programming interface (API) of a TPS (Eclipse, Varian, Palo Alto, CA, USA). A set of beam angle ranking scores (S beam ) for all angles is predicted using the deeplearning model. The quality of the treatment plans created with the predicted S beam guidance was statistically compared for equi-spaced and clinical plans. The evaluation was performed using the dose-volume histogram (DVH) parameters.

Patient Database
Patient data from 50 liver cases, consisting of average intensity projection (AIP) CT images calculated from the fourdimensional (4D) CT images of 40-60% phases, a digital image communication in medicine-radiation therapy (DICOM-RT) structure file, and a DICOM-RT plan file, were used in this study. The patients were originally treated using the proton DS delivery technique configured with three DS fields at the National Cancer Centre in the Republic of Korea (40).
To automatically access information of interesting structures in the TPS, an in-house software was developed using the Eclipse script API. The eight structures of interest included the body, total liver volume (TLV), primary gross tumor volume (PGTV), duodenum, stomach, esophagus, heart, and spleen. The body contour included areas such as immobilizers that should be considered for dose calculation.

Geometric Information and a Set of Beam Angle Ranking Scores for Beam Angle Optimization: Input and Output of the Deep-Learning Model
To train a beam angle optimization network (BAODS-Net), geometric information extracted from the AIP CT images and DICOM-RT structures (an input of BAODS-Net) and S beam generated from RT-plan (an output of BAODS-Net) were used.
The geometric information was collected by the ray tracing method (4). A ray, which is a collector, was determined to penetrate from the body contour to the isocenter, and the path of the ray was tracked in a 3D treatment room coordinate system. The ray collected geometric information by dividing the length of 1,000 mm into 4,000 bins ( Figure 1A). The collected data included the following: specifically HU from the AIP CT images as a double data type and the anatomy position information of interest of eight structures as a binary data type ( Figure 1B). Thus, the shape of geometric information extracted from a collector was 9 × 4,000. Geometric information collected by the penetrating ray is useful data to evaluate the best DS field that considers range uncertainties and OAR positions. In this regard, additional 24 parallel rays penetrating the target volume were created. In detail, the positions of the 24 rays were automatically determined at an isometric angle on two ellipses with different radii placed on the PGTV cross-section in the isocenter plane observed from the beam's eye view ( Figure 1C). Thus, at a specific gantry angle, 25 data collectors extracted nine geometric information points ( Figure 1D). In the same manner, the patient CT anatomy and HU data for all directions were collected at angles from 0°to 358°in steps of 2°. Finally, the shape of the geometric information extracted by the 25 data collectors on coplanar was 40,500 × 4,000. The geometric information was reshaped to input 4D tensor (batch, channel, height, and width). For pre-processing, only the HU values were normalized by using Z-score normalization per patient. The Zscore normalization is a standardization method for the fast convergence of deep-learning models. For the case of the anatomy position information of eight structures, normalization was not performed.
A S beam ground-truth for each patient was generated using the gantry angle information in the DICOM-RT plan files with steps of 1°. The S beam comprised one-dimensional (1D) data continuously ranging from 0.0 to 1.0. The gantry angles used in the clinic were assigned a value of 1.0; otherwise, a value of 0.0 was assigned. Then, to induce effective optimization of BAODS-Net, a normalized Gaussian filter was applied. Finally, the size of S beam was 360, and the shape of S beam was reshaped to a batchconsidered shape. An example of a reference S beam is shown in Figure 2.
K-fold cross-validation (CV) could provide a better indication of how well the BAODS-Net was universalized to unobserved data. We performed a patient-wise K-fold CV method (K = 5), and all datasets were divided into five disjointed and identically sized subsets (41).

Double-Scattering Beam Angle Optimization Network
In this paper, the BAODS-Net based on a CNN is proposed as shown in Figure 3. It consists of two main stages: a feature extractor and a predictor. The feature extractor was configured with eight convolution layers to extract distinguishable features of the geometric information by applying a convolution layer with various strides. The first convolution layer was designed with dimensions of 1 × 9 for the width and height, respectively. The convolution layer was operated with a stride size of nine and with the same padding option. This is because the nine geometric features extracted by a ray are intended to be integrated into a weighted geometric feature. The second convolution layer was designed with dimensions of 1 × 25 for the width and height, respectively. A stride size of 25 was used to integrate each weighted geometric feature extracted from the 25 rays into one weighted ray representing a specific angle. The remaining part of the feature extractor was designed with six convolution layers (width: 3, height: 3, stride: 1) and max-pooling layers. The output of the feature extractor was flattened and then passed to the predictor. The predictor was designed with four fully connected (FC) layers to continuously predict S beam from 0°to 359°. Although the FC layer is computationally expensive, it can effectively predict S beam with a non-linear activation function because the FC layer has a structure agnostic property.
Batch normalization was applied for fast convergence in the optimization process and to ensure the robustness of the performance (42). Randomly initialized biases were added for all layers. Additionally, the activation function was the leaky rectified linear unit (43), which was used to maintain the contribution of negative data, and the adaptive momentum estimation optimizer was employed (44). The BAODS-Net trained for 5,000 epochs with a learning rate of 0.001 and weight decay of 0.0002, and to compute the running average of the gradient, b 1 was 0.9 and b 2 was 0.999. The output shape and parameters of the layers composing the BAODS-Net are summarized in Table 1.

Training and Validation of BAODS-Net
In the training process, BAODS-Net was optimized to predict S beam by using the training data. To find the training loss function that could achieve the best BAODS-Net performance, we compared the performances of the three types of loss functions: L1 loss (Eq. 1), L2 loss (Eq. 2), and smooth-L1 loss (Eq. 3). The prediction accuracy for S beam was calculated by using mean squared error (MSE) between the reference and predicted S beam .
L1 loss(x, y) = 1 L2 loss(x, y) = 1 Smooth L1 loss(x, y) where (x, y) is the reference S beam and predicted S beam , respectively, and n is the number of samples. The hyperparameter beta (b) in Eq. 3 was a value for applying additional weight to the loss. b was empirically determined to be 0.5.
In this experiment, the data of the first fold were used. Specifically, the smooth-L1 loss could be interpreted as a combination of L1 and L2 losses.
To improve the BAODS-Net performance and reduce its generalization error, a 1D augmentation technique, which is a 1D data translation ranging from −2 to 2°, was applied to the reference S beam . The augmentation data were randomly generated for each epoch. Through the K-fold CV principle, we independently conducted five different runs for five separate CV datasets to evaluate the BAODS-Net performance.

Plan Creation With BAODS-Net
The procedure for creating a three-ports proton DS plan with the guidance of BAODS-Net (Plan BAO ) was as follows: (i) the patient CT anatomy and HU data were automatically extracted by using an in-house software based on Eclipse API, (ii) the patient CT anatomy and HU data were fed into the BAODS-Net, and then the BAODS-Net output (S beam ) was predicted; (iii) the specific angles in S beam were selected according to a selection rule. The rule preferentially selected the three gantry angles corresponding to the highest score. However, if the interval was less than 30°, the next priority angle was selected; (iv) the collimator and compensator were designed using the default TPS option without manual modification for objective evaluation of BAODS-Net performance, and (v) the field weight was set to one for all fields.

Plan Comparison of BAODS-Net, Equi-Spaced Angle, and Clinical Plan
To validate the liver treatment plan quality for 50 patients, results were obtained by combining the results of the five folds. The DVH parameters were analyzed for Plan BAO , the equi-spaced plan [gantry angles were fixed at 0°, 120°, 240°(Plan Equi )], and the clinical plan (Plan Clinic ). Although the equi-spaced plan is rarely applied in clinics for proton beam by a planner, we added to the equi-spaced plan for comparative study (4). The Plan Clinic s were created by a qualified planner with 5 years of clinical experience. The evaluation metric is defined below, and the conformity index (CI) was calculated for PGTV (Eq. 4).

Conformity index
where TV is the target volume and PIV is the prescribed isodose volume. The closer PIV is to TV, the closer the CI is to 1. The volumes V 27Gy and V 30Gy for TLV as well as the mean, minimum, and maximum doses for OARs were calculated. V xGy represents the volume percentage of the whole organ receiving a dose ≥xGy. For statistical analysis, these results were compared with the paired-samples t-test. All statistical analyses were implemented using SAS 9.4 software (SAS Institute Inc., Cary, NC, USA), and the statistically significant level was set at P = .05.

RESULT Performance Comparison of L1, L2, and Smooth-L1 Loss Functions
To determine the optimal loss function, BAODS-Net was trained using L1, L2, and smooth-L1 loss. The training time was approximately 120 h by using the training data of the first fold, which corresponded to 5,000 epochs when using an NVIDIA Quadro GV100 graphics processing unit (GPU) (NVIDIA, Santa Clara, CA, USA). The seed number was fixed in the training procedure. The MSE between the predicted and reference S beam was evaluated when L1, L2, and smooth-L1 loss were used for each model training procedure. At 5,000 epochs, the MSEs were 0.031, 0.011, and 0.004 for L1, L2, and smooth-L1 loss, respectively. As a result, the smooth-L1 loss was adopted as a metric for the training loss.

Plan Comparison of BAOBS-Net, Equi-Spaced Angle, and Planner
To evaluate the cases of 50 patients as test data, BAODS-Net was trained and evaluated with five different folds, and the training losses in the five different runs were recorded. At the 5,000th epoch, the mean and standard deviation of MSEs between the predicted and reference S beam for the five folds were 0.0037 and 0.0006, respectively.
The plans created by guidance with BAODS-Net, equi-spaced angle, and the planner method were compared using the DVH parameters. In Figure 4, the paired two test cases of Plan BAO and Plan Clinic are visually analyzed, including the reference, predicted S beam , and 2D dose distribution at the isocenter plane. The star markers in Figures 4A, D are the gantry angles finally selected for Plan BAO .
Data from a total of 37 out of 50 patients were used for comparative evaluation; 13 cases were outside the approved proton range of the TPS for using the equi-spaced angle. The errors occurred in the right posterior oblique (RPO) field of 120°.
For the mean dose of OARs, the results of statistical comparison between Plan BAO and Plan Clinic are summarized in Table 2. The results of statistical comparison between Plan BAO and Plan Equi are summarized in Table 3. As a statistical result, the mean dose has no significant differences between Plan BAO and Plan Clinic (P >.05), while the mean dose of Plan Equi has a significant difference with Plan BAO (P <.05). These results signified that Plan BAO is superior to Plan Equi and similar to Plan Clinic in OARs. The mean dose is visualized for each structure in Figure 5 as a boxplot. The central mark (red) indicates the median, and the top and bottom edges of the box indicate the 25 th and 75 th percentiles, respectively. The whiskers (-) extend to the most extreme data points, while not considering outliers (+). Table 4 summarizes the average of the mean, minimum, maximum doses of OARs for the three planning methods. As a result, guidance using the BAODS-Net method may engender a plan with a quality similar to that created by the planner. In the case of the equi-spaced plan, the quality is relatively low compared to that of the clinical plan.

DISCUSSION
The conventional procedure for creating a proton DS treatment plan is time-consuming and planner dependent. BAO can be utilized as a logical step for the development of efficient and optimal proton plans, similar to studies finding optimal fields in the static IMRT planning area. To date, there is no clinically applicable commercial software for BAO or for enabling intuitive comprehension by a planner.
In this study, we designed BAODS-Net, a new deep-learning based method of BAO for proton therapy. BAODS-Net is based on a CNN and employs the patient anatomy and HU data from the DICOM-RT structure file and AIP CT images, which are automatically extracted using the Eclipse API. The output is the predicted S beam as angle ranking information, which could be used as a priori knowledge to guide the determination of the three gantry angles used for clinical practice.
According to the study results, the proposed method produced clinically acceptable and practical plans. The time required to create a proton DS treatment plan was decreased to approximately 5 min; specifically, the predicted S beam was calculated within approximately 0.2 s through BAODS-Net.
This study provides key contributions and is distinguished from recent BAO research in several ways. This study is the first to employ the BAO method of the proton DS delivery technique, and a deep-learning based one-stop solution was developed. By leveraging this solution, a planner can refer to the predicted S beam in the commercial TPS using Eclipse API. In contrast, conventional BAO research required multiple steps to solve the BAO problem, specifically optimizing the fluence map and then computing the dose influence matrices. However, in our research, only the patient CT anatomy and HU data were used for BAO procedures without dosimetric information from candidate beams. Similarly, Barkousaraie et al. (34) proposed a BAO method using the art column generation (CG) method and a CNN. The architecture of their deep-learning model is based on U-net (45), and the model is trained to mimic the result of the CG algorithm by using only extracted features from the patient anatomy. Although this approach also does not directly use dosimetric information for BAO, more time and additional effort are required to obtain the CG algorithm results.
The following factors may have affected the measurement accuracy. The predicted beam could not provide an optimal  beam configuration because a reference S beam was originated from Plan Clinic . In addition, the Plan BAO was generated without manual modification/optimization procedure such as beam weight, collimator design, compensator design, etc. In other words, it means that the plan quality of Plan BAO has the scope for improvement. In this study, we considered only coplanar proton DS plans for liver cases. However, if the search space is expanded for a non-coplanar proton DS plan, the proposed method could be applied to a non-coplanar proton DS plan for other diseases. Meanwhile, it should be noted that the field design of Plan Equi , specifically anteriorposterior, RPO, and left posterior oblique can be disadvantageous for liver cases. However, according to Figure 5 and Tables 2-4, it can be confirmed that Plan BAO can create a plan of similar quality to that of Plan Clinic by considering the patient CT anatomy and HU data.

CONCLUSION
In this paper, we validated the feasibility of using BAODS-Net for BAO of the three-port proton DS plan. BAODS-Net only used geometric information automatically extracted through the Eclipse API and could successfully predict the S beam for the planning. The results clearly showed its potential for facilitating the three-port proton DS planning. The BAODS-Net dramatically reduced the planning time and brought us one step closer to real-time adaptive proton radiotherapy. Finally, the quality of Plan BAO was statistically verified to be similar to that of Plan Clinic in the mean dose of OARs (P >.05)

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: 10.6084/m9. figshare.14034143.

AUTHOR CONTRIBUTIONS
WC and HK contributed to the conception and design of the study. WC, SA, and SJ organized the database. WC, SY, and SL performed the statistical analysis. HK, SM, and TK provided guidance on methodology and the overall project. SBL, DS, YL, and JJ provided lab and technical support. WC wrote the first draft of the manuscript. All authors contributed to the article and approved the submitted version.  A higher score from the comparison of the three methods is highlighted in bold for organs-at-risk. Bold values stand for having statistical significance. TLV, total liver volume; PGTV, primary gross tumor volume.