A Force Feedback Pedical Screw Implantation Surgical Trainer with Automatic Vertebra Segmentation

The spine is one of the most important parts of the human body. Spine diseases owning to aging or injuries would be treated by spine surgeries. However, the surgeries are highly dangerous and the surgical training can increase the success rate. In this paper, an interactive surgical simulation system for the training of pedicle screw implantation is introduced. We propose an automatic and efficient method to segment the vertebra from the CT images for building the vertebra model. The mean shift method which can cluster points having similar property into groups is used in data clustering. The clustered groups related to the bone region are chosen as the initial model. And then, we apply the Chan-Vese level set method which refines the vertebra contour by using the concepts of statistics. After segmentation, we then apply the marching cube algorithm to generate the 3D vertebra model, and pass it to a virtual surgical environment. A haptic device is adopted to interact with this model to fulfill the surgical simulation in the virtual environment. The average dice coefficient of the segmentation results is 0.92.


Introduction
The spine is one of the most important parts in the human body, it plays a significant role in the human daily +This work was supported by MOST, Taiwan under grant MOST-103-2221-E-006-027-MY2 activities such as walking, lifting objects, or even just standing.The main structure of the spine is composed of a series of separated vertebras and each pair of neighboring vertebras are connected with an intervertebral disc which acts as the function of cushion.But the intervertebral disc will be worn or torn during the aging of human beings.As a result, the vertebras may collide without buffering and lead to spine diseases.Spine repair surgeries therefore become necessary in some of the diseases.However, the spine surgeries are highly critical and the surgical training can help surgeons to elevate the accuracy and success rate of surgeries.In this paper, an interactive simulation system pedicle screw implantation is introduced.
In order to accomplish the surgical simulation system, we have to construct the model based on the vertebra segmentation result from medical images and build a virtual surgical environment.As the preprocessing step, the CT images are smoothed by an anisotropic diffusion (AD) filter which can remove noises while preserve the significant parts of the image.Next the mean shift clustering is applied to create automatically an initial model which is a good prior to the segmentation and can help the refinement step to converge quickly.Finally the Chan-Vese level set method is applied to refine the contour.After segmentation, we apply the marching cube algorithm to generate the 3D vertebra model and construct a virtual environment.A haptic device can interact with this model in the virtual environment to accomplish the surgical simulation.And the on-line force feedback is calculated during the interaction.
The segmentation of vertebra from the CT data images is an important step to build the model for the surgical simulation system.There are many studies to do segmentation by different image processing methods.Chien (1) generated an initial vertebra model by thresholding the CT images and a finer reclassification was used to improve the segmented boundary.However, a blurry boundary region may affect the performance of segmentation.Naegel (2) applied a watershed method to do segmentation with markers, and various steps were needed to find the markers when the image contents were complex.Freedman (3) applied a graph cut method in which a level set shape prior was needed.Chen (4) applied a level set method to do segmentation where a targeted initial contour was needed.Lim (5) applied a 3D level set method with Willmore flow as a prior but the computational cost was very high.Hsiang (6) applied a 3D level-set method with an prior initial model generated by thresholding which was not sensitive to both clear and blur regions.Surgical simulations are very helpful on increasing the success rate, Manbachi (7) made a review on pedicle screw techniques and trainings.Pan (8) introduced a simulation of laparoscopic rectum surgery with force feedback device.Tsai (9) provided a simulation of spine surgery with force feedback device in which cutting simulation was mainly concerned.

Data Set
We use the computed tomography (CT) images of the spine as the data source.There are three L5 vertebra data sets with axial, sagittal and cornel views.The axial view data are used in the segmentation, and the resolution of each image is 512x512 pixels, each data set contains about 100~300 images.The pixel spacing is 0.50mm x 0.50mm and the slice thickness is 0.6mm.Each segmentation result also includes the neighbor vertebras of the target L5 vertebra.

Preprocessing
It is because the CT raw images contain significant noises, so we need to remove the noises by smoothing the image in order to get better segmentation result.The anisotropic diffusion (AD) filter which is also called the Perona-Malik diffusion filter (10) is applied to remove the noise on CT image with good preservation of image features like edges, lines.For our data sets, we apply the 2D AD filter on each image.Figure 1 shows the original image and the image after AD filtering.

Initial Model Construction
After the preprocessing step, we get the smoothed and edge-preserved images.Our target here is to perform an adaptive algorithm to roughly separate the foreground from the image in a fast manner where the vertebra is regarded as the foreground.And then the foreground result is then passed to the refinement segmentation step which is the Chan-Vese level set process.In order to automatically generate this contour for the final segmentation step, the mean shift clustering algorithm (11,12) is applied here.
The mean shift clustering is used to partition data into groups where each group contains a data center with a peak value of kernel density estimation.The kernel density estimation is given by Where n is the number of data, h is the window size of the kernel K, where h d is the d-th power of h, and x is the data at the center of the kernel and ϵ .In this study, the Epanechnikov kernel (13) is applied.
Because the mean shift clustering is an iterative process to find the cluster center, we first give an arbitrary point as the cluster center, at every iteration it moves towards the higher density place.The final cluster center locates at the position where 0. So we can try to observe the change of to find the convergence point of the estimation.The gradient of the kernel density estimation is given by is the normalization constant.For the Epanechnikov kernel, k is the profile of K and is given by The equation ( 2) can be expressed into two parts.The first part is a scalar, and the second part is the mean shift vector and is given by In our case, the data dimension d=3, and the top left pixel in the image becomes the first arbitrary center.The mean shift vector is calculated by equation ( 4) and the data inside the kernel are labeled as a group belongs to the same cluster.It is because the mean shift vector always points to the direction of maximum increase of the kernel density estimation.This mean shift vector is added to the current center to form a new center.The shifting process is repeated until the length of the vector becomes stable or the maximum number of iterations is reached.After a cluster is found, a new arbitrary center is assigned sequentially on the image and the shifting process repeats.The mean shift clustering finishes when every data point is clustered.
The mean and standard deviation of the intensity of the image are calculated based on the resulting clusters.Finally, the clusters which have the average intensity higher than the mean plus 1.75 times of the standard deviation will become  And the initial contour generated by the Otsu's thresholding contains some undesired tiny regions is shown in Figure 2(d) for comparison.

Image Segmentation
The initial model generated is not smooth and serrated as shown in Figure 2(c), so this is refined by the Chan-Vese (CV) (14) level set method.The CV level set method can find the contours which is hard to detect by the gradient method.The CT images of vertebra include a number of slices in the axial view.If the slice is taken on the middle of the vertebra body, the contour of the vertebra is clear as in Figure 3(a).However, if the slice is taken near to the top or bottom of a vertebra body, an ambiguous region and some blurry edges are observed as in Figure 3(b).The traditional level set method is usually constrained on the blurry cases.In order to conquer this, we apply the CV level set method.
The CV level set method embeds the concept of the statistics concerned with the homogeneity of the foreground and background.The level set function can be initiated in many ways.In this study, we generate the initial model based on the distance transform on the obtained contour in the previous step.This transformed model meets the level set definition: the value of the foreground is higher than zero, the value of the background is smaller than zero and the value of the contour is zero.The energy function of CV level set is given by Where C is contour, , is the average intensity of the foreground and background respectively.Length and Area are the contour length and area of the contour respectively.f(x) is the original image., , , are four weighting parameters.The evolution is performed by the differential of the energy.The best contour can be found by minimizing this energy.In our study, the evolution stops when the level set map becomes stable or the maximum iteration number is reached.Figure 4 shows the initial contour and the result of CV level set where the contours are shown in the red color.

Environment Setup
Our surgical simulation system is built by Visual C++ with OpenGL 2.7, and the Geomagic Touch (Sensable Phantom Omni) (15) is used as the force feedback device.The device has 6 degree-of-freedom and three directions of force feedback.The surface model is constructed on the segmented vertebra volume by using the marching cube algorithm (16) .Because of the number of vertices and faces generated by the original marching cube is too large, it will cause too much computation when rendering.The HC Laplacian (17) filter is applied to smooth the surface model.We thus place a simplified model with similar appearance to increase the efficiency of system interaction in the virtual environment.The end point of the device in the real world is displayed as an interactive point in the virtual world with an update rate of P Hz.The tool model transformation is to multiply the transform matrices with the original tool model.Figure 5 shows a scene of the tool interacting with the vertebra model.The force applied real-time on the vertebra in the virtual environment is fed back to the haptic device during the tool interacting.

Force Simulation
The intensity of the original image is used as the reference information of force feedback.The vertebra is constructed by the spongy bone and the cortical bone.The cortical bone is denser and has a higher intensity value in the CT image, and generates a greater force feedback.In our system, we combine the surface model and the original image data together to accomplish this idea.
When calculating the force on the device, we extend a solid sphere from the interactive point, the voxels of this sphere take part in the calculation of feedback force.We build the force model (18) as shown in Figure 6 in which the circle is the cross section of the sphere, is the tool vector, the yellow colored voxels are the surface voxels of the sphere and the green colored voxels are the voxels inside the sphere.The calculation of the force of each voxel is given by * * * , Where k is the intensity of each voxel, D is the ratio of the numbers of bone to non-bone voxels of the sphere.is different with respect to the voxels at the sphere surface or inside the sphere.For the voxels inside, the force vector points to the center of sphere.For the voxels on the sphere surface, the force vector points to the center of the tool. is the angle between the tool vector, which is the axis of the tool, and the force generated by each voxel.Finally, the force feedback is calculated by summing these forces generated by the voxels of the sphere.

Results and Discussion
For the segmentation results, Figure 7(a), (c) and (e) show the results where the slices are taken at the middle of the vertebra, and Figure 7(b), (d) and (f) show the results which are difficult to accurately obtain when the slices are taken between the intervertebral disc and the end of a vertebra.It can been seen that our segmentation results are good on both normal and difficult cases.
The mean shift clustering method is applied to generate the initial contour for Chan-Vese level set segmentation in this study.As the initial contour can also be generated by different methods, the initial contours generated using a circle placed at the image center or by using the Otsu's thresholding method are compared with the proposed method for their level set segmentation results as shown in Figure 8 and generated by the Otsu's thresholding and our proposed method are both close to the ground truth but is more fitting by our proposed method.The segmentation results in the right column for the initial contour generated by the Otsu's thresholding method is too much over-segmented and our proposed method is better fitted.However, there are still two small unwanted regions and is a little bit over-segmented in the resulting segmentation of the proposed method.
The Dice coefficient (DC) is a common and easy way to examine the accuracy of segmentation.We have 3 data sets, our evaluation is to randomly choose 5 images from each data set and calculate the DC from our results with the ground truth which was drawn by an expert.The DC formulation is given by Where A is the area of the ground truth, B is the area of the segmentation result.The segmentation results are built as 3D models by using the marching cube algorithm.Figure 10  not only a virtual surgical scene, but also the force feedback to real world and the scenario according to the expert's advisement such as the surgery process are designed.And we are expecting the user can be trained by this system to experience the real surgical feeling and give us feedback to improve the system.Figure 11 shows a user using the system.

Conclusion
We build a new force feedback surgical simulation system for the pedicle screw implantation surgery.An automatic and efficient segmentation method has been proposed to segment vertebra and build the model for the system.It does not need to set the initial model point manually and the initial contour is automatically created.After getting the contour which is close to the final contour, we apply the CV level set method to iteratively refine the contour with the statistic concept of homogeneity.Then the vertebra model is constructed with the segmented result by using the marching cube algorithm.A haptic device is then adopted to interact with this model to simulation the surgery.The force feedback is calculated by using the model with the original intensity data to increase the reality.By referring to the experts' opinions, the training system is designed with the scenario of surgical reality.In the future work, we will invite more users to adopt more feedbacks for improving the force interaction and system performance.

Fig. 2 .
The mean shift clustering.(a) Input image.(b) Clustered image.(c) Initial model chosen by mean shift.(d) Initial model generated by Otsu's thresholding.

Fig. 3 .
Images taken at different positions of vertebra.(a) Clear vertebra contour.(b) Blurry vertebra contour.the initial model.Figure 2(a) shows the input image which is AD filtered and Figure 2(b) shows the mean shift clustering result.The initial contour chosen is shown in Figure 2(c).

Fig. 5 .
Fig. 5.The scene of a tool interacting with vertebra.

Figure 9 .Fig. 7 .
The segmentation results with the initial contours a circle placed at the image center show too much over-segmentation on both cases.The segmentation results in the left column for the initial contour The segmentation results (in red color) for different cases.
(a) shows a L5 vertebra model at the front view and Figure 10(b) shows the same model at the back view, the neighbor vertebras are shown together here.We use the segmentation result to generate a surface model and apply it in the surgical simulation.In this system, (a) (b) Fig. 10.The 3D vertebra models.(a) Front.(b) Back.

Fig. 11 .
Fig. 11.The user interface of the system.

Table 1
shows the average results of each data set.Data 2 is an aged patient who has low bone density which causes low intensity in images and affects the segmentation results.Fig.8.Initial contours (in red color) by different methods.Fig.9. Binary segmentation results by different initial contours.(Red contour is the ground truth.)

Table 1 .
The segmentation results of the data sets.