# Feature-Oriented CBCT Self-Calibration Parameter Estimator for Arbitrary Trajectories: FORCAST-EST
## 2 Materials and Methods
### Algorithm  
An overview of the estimation process is in Figure 2, and this will be explained in more detail in the following section. The algorithm is initiated with the acquired images, a prior CT, and geometry information about the size of the detector, as well as the distance to the source and to the isocenter. The first step of the algorithm is to generate simulated projections in a regular grid around the center of the prior CT image. The grid has 95 points each for the rotations around the x- and y-axis, which results in a grid with a spacing of 3.8\({}^{\circ}\). Further, three different detector rotations are used, 0\({}^{\circ}\), 120\({}^{\circ}\), and 240\({}^{\circ}\). On each of these images, the AKAZE algorithm detects features and extracts the feature descriptors. The simulated projections are then deleted, and only the feature points, descriptors, and projection rotations (\(\alpha\), \(\beta\), and \(\gamma\)) are saved. They are calculated once and then used for all further calibration. To find the approximate position of an acquired image, the algorithm first detects the features. They are then matched with each set of features from the simulated grid  
Figure 2: Pipeline of the estimation algorithm. Rectangles are the processes, and the trapezoids are the store data. projections and the matched feature points are counted. To save time, the algorithm operates along this grid with a step size of four, and it then selects the five grid points with the most matched feature points. Next, the grid points surrounding these five points are compared to the acquired image in the same way. For the grid point that has the most matched pairs, the projection rotations are returned. These are the current approximation of the rotations \(\hat{\alpha}_{1}\), \(\hat{\beta}_{1}\), and \(\hat{\gamma}_{1}\) (Figure 3a). This selected grid position is most likely the closest to the target position, but it can also be a projection from the opposite side. These projections, simulated from the wrong side, will be corrected later. Before that, the detector rotation is approximated. The average value of the feature point coordinates is used as the center point; this is performed separately for the real and simulated image. Then, the coordinates of the feature points are converted to polar coordinates by using the averaged center point as the zero point. The angle for each feature point in the simulated image is subtracted from the angle of the matching feature point in the real image. The median of these differences is the new approximated detector rotation \(\hat{\gamma}_{2}\). Listing 1 shows this in pseudocode. ```
1 defapproximate_detector_rotation(current_parameters):
2#simulateprojectionandtrackfeatures
3 simulated_projection=ForwardProjection(current_parameters)
4 simulated_feature_points=trackFeatures(simulated_projection,
5 real_projection)
6#calculatecenterpoint
7 simulated_mid=mean(real_feature_points,axis=0)
8 real_mid=mean(simulated_feature_points,axis=0)
9 sim_points=simulated_feature_points-simulated_mid
10 real_points=real_feature_points-real_mid
11#calculateangle
12 angles=(arcand2(sim_points[:,0],sim_points[:,1]))
13-arctan2(real_points[:,0],real_points[:,1]))*180.0/PI
14 angles[angle<-180]+=360
15 angles[angle>180]-=360
16 detector_angle=median(angles)
17#testinwhichdirectiontotrotate
18 proj=ForwardProjection(applyRotation(current_parameters,
100,-detector_angle))
19 points=trackFeatures(proj,real_projection)
20 diffn=points-real_feature_points
21 proj=ForwardProjection(applyRotation(current_parameters,
2200,+detector_angle))
22 points=trackFeatures(proj,real_projection)
23 diffp=points-real_feature_points
24 ifsum(abs(diffn))<sum(abs(diffp)):
25 return-detector_angle
26 else:
27 returndetector_angle
```  
**Listing 1**Approximation function for the detector rotation. A similar approach is taken to correct the projections taken from the opposite direction.