Improving pedicle screw path planning by vertebral posture estimation

Objective. Robot-assisted pedicle screw placement in spinal surgery can reduce the complications associated with the screw placement and reduce the hospital return counts due to malfunctions. However, it requires accurate planning for a high-quality procedure. The state-of-the-art technologies reported in the literature either ignore the anatomical variations across vertebrae or require substantial human interactions. We present an improved approach that achieves pedicle screw path planning through multiple projections of a numerically re-oriented vertebra with the estimated posture. Approach. We proposed an improved YOLO-type neural network model (YOLOPOSE3D) to estimate the posture of a vertebra before pedicle path planning. In YOLOPOSE3D, the vertebral posture is given as a rotation quaternion and 3D location coordinates by optimizing the intersection over union of the vertebra with the predicted posture and the actual posture. Then, a new local coordinate system is established for the vertebra based on the estimated posture. Finally, the optimal pedicle screw path trajectory is determined from the multiple projections of the vertebra in the local coordinates. Main results. The experimental results in difficult cases of scoliosis showed that the new YOLOPOSE3D network could accurately detect the location and posture of the vertebra with average translation and orientation errors as small as 1.55 mm and 2.55°. The screw path planning achieved 83.1% success rate without breaking the pedicle cortex for the lumbar vertebral L1–L5, which is better than that of a doctor’s manual planning, 82.4%. With the clinical class A requirement to allow less than 2 mm out of the pedicle cortex, the success rate achieved nearly 100%. Significance. The proposed YOLOPOSED3D method can accurately determine the vertebral postures. With the improved posture prior, better clinical outcomes can be achieved for pedicle screw placement in spine internal fixation procedures.


Introduction
Spinal degeneration is a common disease in the aging population, while spinal fractures and deformities may involve all age groups due to trauma and tumors, etc (Chin et al 2007, Tomé-Bermejo et al 2017, Zou et al 2020. A surgical procedure with pedicle screw fixation is one of the commonly used solutions for the treatment of these diseases (Cordemans et al 2017, Kong et al 2017, Zhang et al 2017, Hou et al 2018. Improper screw placement can cause serious clinical complications as the incidence rate of non-guided surgery can reach as much as 52.58% (Battié et al 2012, Lau et al 2013, Aoude et al 2015, such as neurovascular injury, dural sac injury, organ and muscle injury, paralysis, blood loss and even death. The use of assistive techniques with navigation to place pedicle screws has been a major trend in the practice of spinal surgery (Driver and Groff 2021). Compared with freehand procedures, the navigated techniques can improve the accuracy of pedicle screw placement, thereby reducing the complications associated with screw placement and reducing the number of returns to the operating room due to malposition (Overley et al 2017, Du et al 2018, Vaishnav et al 2019, Fatima et al 2021, Roser and Hebela 2021, Tarawneh and Salem 2021. Several studies have shown that the accuracy of pedicle screw fixation assisted by a robot system is better than that of manual surgery guided by fluoroscopy, with an acceptable success rate of 94%-98% (Roser et al 2013. Moreover, as a result of training, inexperienced doctors are more likely to perform the procedure with higher accuracy and less operation time (Siddiqui et al 2019, Morse et al 2021, Yu et al 2021. However, the current commercial robotic surgery systems usually lack automation as the doctors must manually plan the direction and position of the screw track on the preoperative or intraoperative images before operation (Kochanski et al 2019, Scalia et al 2020. Current pre-surgical path planning uses a manual drag-and-drop strategy on an intraoperative cone-beam CT (CBCT) image displayed on a visualization workstation to determine the screw trajectories with repeated modifications until a satisfactory result is reached. This process will interrupt the workflow and may cause human errors. Therefore, various computer-aided planning methods have been introduced intending to improve surgical outcomes. Wicker and Tedla (2004) and Lee et al (2011aLee et al ( , 2011bLee et al ( , 2012 planned the optimal screw trajectory by manually extracting the pedicle region and reconstructing the geometric shape and characteristics of the pedicle. Knez et al (2019), Mischler et al (2020), Wi et al (2020), Soliman et al (2021) and Caprara et al (2022) calculated the optimal screw trajectory using bone mechanics or bone mineral density characteristics obtained from CT images. Daemi et al (2015) and Solitro and Amirouche (2016) calculated the preoperative screw trajectory based on the correlation between the screw placement trajectories at different angles and vertebral bone mineral density. Knez et al (2016), Xiaozhao et al (2016) and Knez et al (2018) developed a method to automatically segment vertebral bodies and pedicles. The trajectory planning was performed based on the geometric features of the pedicle. Goerres et al (2017) and Vijayan et al (2019) used the atlas method to align the screw path and vertebra on the known vertebra to the target vertebra for screw trajectory planning. Cai et al (2019), Kausch et al (2021), Ma et al (2022) and Qi et al (2022) developed the methods of spine segmentation and preoperative screw trajectory planning based on artificial intelligence (AI).
The above works have obtained reasonable experimental results, but there are still some drawbacks. For example, these methods either require the operators to manually select vertebral body information from CT images, or consider the holding force and stability of screws but ignore the geometric characteristics of the vertebra. The methods based on AI and atlas require a large number of atlases of vertebrae containing screw tracks, and their screw planning methods do not consider the differences between different vertebral structures, which may be problematic in a personalized planning method for different patients.
This study presents an optimal pedicle screw trajectory planning algorithm reinforced by vertebral posture estimation, which is suitable for a personalized screw trajectory plan and without manual intervention. In the optimal screw trajectory planning, it is important to mark the avoidance region as a limiting force in costfunction-based mathematical modeling. The procedure of marking the avoidance region can be completed with a contour recognition (Suzuki et al 1985) with the assistance of binary segmentation of the vertebra. The location and posture of the segmented vertebra are estimated by a convolutional neural network model from preoperative CT volume data. Then, a local coordinate system can be established for individual vertebra, and an optimal screw path can be estimated from multiple directions based on the local coordinate system.

Materials and methods
An overview of the method for automatic planning of the screw trajectory reinforced by posture estimation is shown in figure 1(a). The location and posture of a vertebra can be estimated by a convolutional neural network model from preoperative CT volume data. Next, a local coordinate system can be established for that individual vertebra. Then, the screw trajectory can be determined based on the geometric characteristics of the vertebra in the local coordinate system. In implementation, a segmented binary vertebra image from the U-NET model is used in trajectory planning for better recognition of the geometric features of the vertebral body.

Convolutional neural network architecture
In this work, we introduce a modified Yolo network, YOLOPOSE3D, to infer the localization and orientation of the vertebra from CT volume data. The model contains a network backbone and a prediction header like the YOLO v3 model. The backbone extracts the image feature and the header infers the localization and orientation of the target from the extracted feature ( figure 1(b)). In our study, we choose the bonding boxes of (70, 60, 50), (120,135,75) and (130,90,80), with the numbers in parentheses indicating the width, height and length in millimeters as the priors of the vertebral sizes in three dimensions to reduce the computational load.
Following Yolo v3, the YOLOPOSE3D predicts the center coordinates and the dimensions, p x , p y , p z , p w , p h and p l for each bounding box (dashed red box in figure 2) as: where c x , c y and c z are offset values from the top left back corner of the image to the voxel cell, where the target center is located, and t x , t y , t z , t w , t h and t l are the adjustment parameters for the coordinates inferred by the network model. The  w ,  h and  l are the width, height and length of the prior box (dashed dark blue box in figure 2). Figure 1. A schematic chart of the proposed workflow for the screw planning process. The outputs t x , t y , t z , t w , t h , t l and * q 1 , * q 2 , * q 3 , * q 4 , can be transformed into vertebral position, posture, and size using (1), (2), and (5).
Generally, the network model can easily learn the target position in Euclidean space, as illustrated in figure 2. However, the task for determining orientation is more complex for a neural network model. The common methods for orientation representation include Euler angles, axis-angle, rotation matrices and quaternions. However, for the Euler angle, it is a periodic function with the period 2π, i.e. there are multiple values representing the same angle. A similar problem exists for the axis-angle representation. Therefore, both the Euler angle and the axis-angle can cause misrepresentation in network learning due to their non-uniqueness issues. For rotation matrices, over-parametrization in rotation regression is a big challenge for neural networks. In addition, the intrinsic orthogonal property for the rotation matrix is difficult to enforce in neural network training. Inspired by Kendall et al (2015) and Kendall and Cipolla (2017), we use a quaternion as our orientation representation as it can be easily mapped into legitimate rotations by normalizing them to a unit vector shown in equation (2), and it is a continuous and smooth representation of rotation. Generally, in a sparse detection task, the number of negative detections is far greater than the number of positive detections. We desire that the neural network is more sensitive to the positive samples. To achieve this purpose, we use the balanced binary cross entropy function optimization model to predict the target: where the γ p and γ n are the balance factors, and y and y¢ are the ground truth and predicted values, respectively. In this work we set γ p = 4.31 and γ n = 0.01.

Loss for setting improper bonding box parameters
Generally, computing the loss function for a proper bonding box,  box , with arbitrary orientation is a challenge because the overlapping area between two rotated boxes is an undetermined polygon whose area cannot be calculated by a differentiable computing model, such as in a neural network. The conventional solution is to disintegrate the  box into two independent components ( IoUn for a non-rotated box and  rot ) to regress the localization, size and posture for the bounding box, respectively.
where λ box and λ rot > 0 are the scale factor to balance each part of the loss function  box ,  IoUn is the intersection over union (IoU) for two non-rotated boxes, which can be obtained by simple calculation. Inspired by Yang et al (2021, 2022), a more appropriate method to calculate  box is given as follows. According to Yang et al (2021Yang et al ( , 2022, an arbitrary-oriented bounding box can be modeled by Gaussian distribution ( )  , m S , expressed as equation (5): where  is the rotation matrix. The volume of the rotated box can be expressed as: Then, the volume of the rotated box for the ground truth g u and for the predicted t u can be defined by With Gaussian product lemma, the distribution function of the overlapping area, which is proportional to the overlapping area for two rotated boxes if the center points are superposed: Then, the overlapping area Iou u can be defined in the way of equations (6) and (7). It is clear that the overlapping volume is only related to Σ g and Σ t . Therefore, to achieve the final IoU-based loss function of rotated boxes,  IoUr , we introduce the distance factor to solve this problem: We also introduce  rot to regress the posture for the bounding box as  IoUr cannot directly regress the posture for the bounding box. Then, the  box can be expressed as:

Overall loss function
Finally, the overall loss function includes the loss for the vertebra target (i.e. the object loss,  obj ) and the loss for setting improper bonding box parameters (i.e. the bounding box regression loss,  IoU , and the rotation loss,  rot ): where λ obj , λ box and λ rot > 0 are the scale factors to balance each term of the loss function , and  Iou represents  IoUn or  IoUr . We set λ obj = 1, λ box = 1 and λ rot = 1.

Vertebra segmentation
A binary vertebra image is required to identify the geometric features of the vertebral body. Therefore, we segment the vertebra with a 3D U-NET (Ronneberger et al 2015) model. In this procedure, we trained the 3D U-NET in a dataset which contains 95 vertebrae of L1 to L5 from 19 patients with doctor-labels. The dataset is split into training data, testing data and validation data with portions of 80%, 10% and 10%, respectively. The average Dice similarity coefficient (DSC),accuracy, precision, and re-call rates are 0.88, 0.99, 0.96 and 0.91, respectively, in the validation results. Figure 3 indicates the segmentation results, which achieve comparable results to the doctor's labeled results. is the unit vector for the screw. The optimal screw path can be found by:

Estimate screw trajectory from projections
The optimal screw placement trajectory in a projection plane can be expressed as a linear equation, as shown:  b . 0 w + = , which maximizes the distance between the point on the avoidance area of the projection plane and the trajectory, where ω and b are parameters of the line.
Assume P θ is the projection at θ angle for the CT volume; the optimal screw placement trajectory in P θ can be determined by ω and b, where ω and b are parameters of the screw trajectory .
is the boundary of the projected M on P θ . The regions split by M¢ can be marked by -1 or +1. It can be expressed as This process can be performed on the binary projection plane using contour recognition (Suzuki et al 1985), which is seamlessly integrated into computer vision libraries, such as the OpenCV library, for convenient implementation. Employing the contour recognition algorithm, we can extract the contours of the avoidance area on the projection, designating the pink contour as −1 and the orange contour as +1, as depicted in figure 4. Then, the optimal screw placement trajectory in P θ can be found as: The above constrained optimization problem can be converted into an unconstrained optimization by the Lagrange multiplier: where α is the Lagrange multiplier with α 0. It can be solved using a support vector machine approach (Platt 1998). Assume P¢ q is the projection at θ angle for an optimal screw placement trajectory with the points on the optimal screw placement trajectory marked 1 and the others marked 0, which can be estimated by equations (13) and (14). The 3D screw trajectory can be estimated in a CT reconstruction style:

Estimating screw trajectory in simplified procedure
In practical deployment, we can estimate the screw trajectory by only two projections since the screw trajectory is simple enough as a line segment. In this work, the projections in the axial (P axial , θ = 0) and sagittal ( ) The x, y and z represent coordinates of the CT volume.

Datasets and processing settings
We validated the proposed method using three sets of clinical data (CT images), which were obtained from Beijing Jishuitan Hospital. The typical dimensions of the axial image are 512 × 512, with a pixel size of 0.5 mm × 0.5 mm and a slice thickness of 0.8 mm. The pixel size and slice thickness may vary with the field of view and slice settings. Dataset 1 contains 76 lumbar vertebrae data from L1 to L4 without a doctor-labeled optimal screw trajectory (DLSTs). The dataset is augmented by randomly rotating (x, y, z-axes, −30 ∼30°) and translating (x, y, z-axes, −10 ∼10 mm) each vertebra. Finally, 1520 re-oriented vertebrae are simulated in dataset 1. Dataset 2 contains 248 lumbar vertebrae data from L1 to L5 with scoliosis. Dataset 3 contains 150 lumbar vertebrae from L1 to L5 with scoliosis. Dataset 3 has the doctor's pedicle screw trajectories (DPSTs) for each vertebra. DLSTs are labeled by an expert physician (with more than 5 years of experience), and subsequently reviewed and amended by two senior physicians.
In the training and validation procedures, we first trained the YOLOPOSE3D in dataset 1 with training 80%, testing 10% and validation 10%. Then, by migrating the training outcome, we further trained the network on dataset 2 with training 80%, testing 10% and validation 10%. Finally, we evaluate the trained YOLOPOSE3D network with dataset 3. We used the following criteria to evaluate the quality of the screw trajectory: deviation angle (DA), point-to-point distance (PPD) and the amount of screw breach through the pedicle cortex (DSOS). For comparison, we implement the method by Kausch et al (2021), which estimates the mask of the screw trajectory using the U-Net network framework. The screw parameters are derived from the predicted screw mask.
All the computations were performed on a HP workstation with dual nVIDIA RTX 8000 GPU and Intel Xeon Gold 6240 CPU. The Visualization Toolkit (VTK) was used to visualize the results. The customized convolutional neural network model with YOLOPOSE3D inscription was written using Pytorch V1.11.0 +cu113. All other programs were written using Python V3.7. The input volumes were 96 × 384 × 384 with intensity normalized to 0 ∼ 1. A batch size of 8 was chosen, the number of epochs was 600 and an Adam optimizer was used at a learning rate of 1e-4.

Validating YOLOPOSE3D
The framework of YOLOPOSE3D was trained and validated with dataset 1 using both the conventional loss function (CLF),  IoUn , and the improved loss function (KLF),  IoUr . With the CLF, the mean absolute errors (MAEs) for translations along the x, y, and z-axes are (T x , T y , T z ) 0.79 mm, 0.87 mm and 0.70 mm, respectively, while the rotations around the corresponding axes (j x , j y , j z ) are 2.48°, 2.25°and 2.28°, respectively. With the KLF, the MAEs for the translations are 0.53 mm, 0.53 mm and 0.49 mm, respectively, and for the rotations are 2.09°, 1.97°and 2.14°, respectively. More information is shown in table 1.
Therefore, the KLF leads to better precision of position and posture estimations compared with the CLF while reducing the standard deviation of the errors. Figure 5 shows the error distributions for the KLF versus the CLF.
Further, we migrated the trained YOLOPOSE3D to dataset 2 and trained/validated it using similar procedures. The MAEs for the translations with the CLF were 0.78 mm, 1.87 mm and 0.97 mm, respectively, and 2.82°, 2.63°and 2.55°, respectively, for the rotations. With the KLF, the MAEs were 0.70 mm, 1.55 mm,   0.56 mm, 2.55°, 2.14°and 2.17°, respectively. We noted that the validation errors in dataset 2 are higher than those of dataset 1. This is because dataset 1 is simulated data and the vertebrae are relatively normal in shape. Dataset 2 is actual data with scoliosis. As table 2 shows, similar to the test on dataset 1, the KLF can achieve better precision than the CLF in both position and posture estimations for scoliosis data. Figure 6 shows the distributions of the errors for both loss models. Therefore, KLF-based YOLOPOSE3D will be used in the rest of the paper.

Effect of YOLOPOSE3D in screw trajectory estimation
We estimated the screw trajectories on vertebral data used in  with the additional information of vertebral position and posture from YOLOPOSE3D. The results were compared with the trajectory estimation results in , where the vertebral posture information is not used. Figure 7 shows the distribution of the quality measure, DSOS, from the previous method ( figure 7(a)) versus the outcome of the proposed method ( figure 7(b)). It shows that the YOLOPOSE3D can substantially reduce the number of screws encroaching the pedicle cortexes (the samples appear in the blue clinically acceptable encroachment (CAE) region).

Screw trajectory estimations for scoliosis vertebrae
Dataset 3 is the scoliosis vertebrae data with doctor's manually planned screw trajectories. With this dataset, we confirmed the effectiveness of YOLOPOSE3D in the screw trajectory planning process. Figure 8 demonstrates  Table 2. Comparison of the position/pose estimation errors for different loss functions for dataset 2 (mean ± SD). the distribution of the DAs between the corresponding estimated pedicle screw trajectories (CPSTs) and the DPSTs for L1 to L5. The average DAs are 4.77°, 5.53°, 5.19°, 4.11°and 7.61°for L1-L5, respectively. Three qualified doctors had confirmed that our results were more clinically acceptable than those of the original doctor's results. For each pedicle screw, we virtually split its full length into 19 equal length segments with 20 points (including two end-points). We investigated the PPD between the CPSTs and DPSTs for all L1-L5 in figure 9 with figure 9(a) for L1-L4 and figure 9(b) for L5. For L1-L4, the maximum average PPD is 2.47 mm, which is located at the anterior end of the screw. The minimum average PPD is 0.79 mm, which is located at the center vicinity of the trajectory. The PPD's results showed that the overall difference between CPSTs and DPSTs is reasonably small. Combined with the DA, we can see that the CPSTs and DPSTs are highly consistent for L1-L4. However, for L5, the PPDs show larger deviations due to the larger size of the vertebra allowing more tolerable room for trajectory variations and the interference of the iliac crest near L5. However, most of the CPSTs are still in line with the clinical requirements.
Although the values of the DA and PPD can imply whether a screw is in line with clinical requirements or not, they cannot show whether a screw encroaches the pedicle cortex or not. Therefore, we further evaluate the DSOS measure, the distance from the screw surface to the pedicle cortex. According to the clinical requirements, it is still acceptable if a screw breaks and extrudes the pedicle cortex no more than 2 mm. The violin plot in   figure 10 describes the distributions of the DSOSs for both CPSTs (figure 10(b)) and DPSTs (figure 10(a)) in L1 to L5 with the assumption of the 6 mm screw diameter, most recommended for lumbar pedicles. A DSOS smaller than −2 mm (in the orange region) indicates a Frank penetration (FP); a DSOS between of −2 mm and 0 (in the blue region) indicates that the screw encroaches the pedicle cortex, but is still clinically acceptable; a DSOS greater than 0 indicates that the screw is within the limit of the pedicle cortex. Accordingly, we further define three additional measures for the screw placement quality, i.e. the MGM, the clinical goodness measure (CGM) and the CAE. For the MGM, the acceptable screw placement requires the screws within the pedicle surface limit (DSOS > 0 mm). For CAE, the screws are allowed to break and extrude the pedicle cortex, but no more than 2 mm. Logically, the CGM covers both the MGM and CAE. Figure 5 shows only a small proportion of CPSTs encroach the pedicle cortex (but within the 2 mm range); all the others are limited within pedicle surfaces. The DSOS results indicate that all CPSTs meet CGM requirements, and their MGM rates compared to DPSTs are 63.8% vs 56.9%, 79.3% vs 73.9%, 93.1% vs 93.1%, and 98.3% vs 96.6% for L1 to L4, while they are 81.0% vs 89.7% for L5. The CGMs for both CPSTs and DPSTs are 100%. These numbers indicate that the quality of screw trajectories calculated by the proposed method is higher than that by the doctor's manual method for L1-L4, but with more error for L5 due to the aforementioned reasons. Table 3 shows the MGM, CAE and maximum distances for screws extruding the pedicle cortex (MEPDs) of dataset 3. We notice that the MGMs for L1 and L2 are relatively low. This is because their smaller pedicle sizes result in a higher possibility of being penetrated during the screw implantation with a standard lumbar 6 mm screw. However, the proposed approach can maximally prevent the screws from penetrating the pedicle cortex (as shown in figure 11). This is also the reason why the MGMs for CPSTs are better than those of DPSTs for L1-L4. A similar trend can also be observed with the quality measure of the MEPD: its values for CPSTs are considerably smaller than those of DPSTs. Once again, larger errors are observed in L5 due to the variation conditions mentioned previously. However, all the screw trajectories on L5 planned by the proposed method are also within a clinically acceptable range. . We can see that the results from Kausch et al (2021) are the best among the competing studies with zero FPs. For comparison purposes, we reproduce Kausch's method on both the normal dataset (data used in Zhang et al (2021)) and the scoliosis dataset (dataset 3). Table 5 and figure 12 show the comparisons between Kausch's method and the proposed method. For the normal spine data used in Zhang et al (2021), Kausch's method only produced 1 FP result out of 190 screws; however, for the data with scoliosis (dataset 3), 70 FPs were produced out of 300 screws. We repeated Kausch's method with the  vertebrae posture prior (generated by YOLOPOSE3D) in the scoliosis data. It achieved much improved results with only two FP outcomes. For the same dataset, there is no FP outcome with the proposed method. More encouragingly, the number of CAE screws with the proposed method is also substantially less than that of Kausch's method.

Discussion and conclusion
This study presents a vertebral posture estimation reinforced optimal pedicle screw trajectory planning method. The YOLOPOSE3D model proposed in this work can estimate the position and posture of the vertebral body accurately. The estimated position and posture information become the reinforcing prior in the screw trajectory planning procedures. We tested the YOLOPOSE3D on both simulated data and real data. The results showed that YOLOPOSE3D could detect the vertebra accurately and provide tight representation of geometric characteristics for the individual vertebrae. The information obtained from YOLOPOSE3D provided better conditions for determination of the optimal pedicle screw trajectories from multiple projections. The validation tests included screw trajectory planning for both normal and scoliosis vertebral scenarios. The results show that all of the calculated screw trajectories are in line with clinical requirements with substantial higher-quality measurement values (MGM and MEPD) in L1-L4. The values for L5 are small due to its surrounding special structures. However, all the screw trajectories on L5 by the proposed method were within the clinically  acceptable range. The proposed method can make the amount of pedicle cortex penetration as small as possible in the cases in which the pedicle cortex must be penetrated. The quality of the screw trajectories planned in our work is better than that planned in the literature. However, this study also has some limitations. The first limitation is that we do not consider multiple vertebrae joint screw placements in this work. Second, due to calculation costs, the input image size of the YOLOPOSE3D model is relatively small (96 × 384 × 384 in our case), and the input image contains only one vertebra. Third, because the geometric structure of L5 is different from that of L1-L4, its geometric feature extraction is more difficult, and the muscles and iliac crest can interfere with the screw planning process; the quality of the screw trajectories planned in L5 is reduced. We also tried the new approach on available thoracic vertebrae (T10-T12) and received acceptable results. However, since we only have a small number of thoracic data, further investigation is needed for more conclusive evaluations. Cervical vertebrae are much smaller in size compared to the lumbar vertebrae and have more restricted accuracy requirements. The proposed approach is a promising candidate, but a thorough investigation is needed. Lastly, in the actual clinical operation, it is desired to plan the screw trajectory on the intraoperative CBCT image for direct application in surgery. However, due to the low quality of the CBCT, the proposed method is carried out on preoperative CT, which requires a registration step for application in surgery.
In conclusion, this paper proposes an improved pedicle screw trajectory planning scheme that is reinforced by vertebral posture estimation. The experimental results show that the screw trajectories obtained from the new method are within the range of clinical safety requirements and better than the doctor's manual results. The new method is suitable for clinical use, and is a better alternative than current clinical practice.