Keywords

1 Introduction

The cervical spine is a complex system, which consists of vertebrae, intervertebral discs, as well as muscles and ligaments. Cervical spine injuries are a major health concern, as trauma in the neck represents the majority of all spinal lesions [9]. Real-time motion assessments of the cervical spine provide means for understanding the natural neck motion and reveal movement abnormalities associated with spinal injuries or medical conditions, such as whiplash or disc herniation. However, no widely accepted standards exist for diagnosing cervical spine injuries in soft tissue. Such injuries are closely related to mechanical factors and thereby in-vivo studies of the spine’s motion can lead to an improved understanding of these injuries [5, 9]. Digital videofluoroscopy (DVF) analysis techniques, which involve recording multiple X-ray images as video frames, have been demonstrated to successfully investigate the spinal motion [1, 2, 6, 7, 15, 16, 20]. Intervertebral kinematics serve as means to describe and quantify the spinal motion. A commonly used method for calculation of intervertebral kinematics is Distortion-Compensated Roentgen Analysis (DCRA) [8, 11, 14, 16]. DCRA is based on landmarks, which are defined in each of the DVF frames [11]. Definition of landmarks is time-consuming and subjective when performed manually, and small errors in the definition of the landmarks can potentially cause large errors in computation of the intervertebral kinematics [13].

Several semi-automatic methods have been successfully developed to track the lumbar vertebrae [7, 10, 19, 20]. Lumbar vertebrae tracking methods have been based on template matching [7, 10, 20] or active contours [19]. Template matching technique involves defining a template for each vertebra, so that this template is continuously tracked throughout the video frames, using cross-correlation as a similarity measure [3, 7, 10, 12, 20]. These methods consider the vertebrae to be rigid bodies and thereby assume that no out-of-plane coupled motion is present. On the other hand, active contour methods do not assume vertebral rigidity [8]. Wong et al. [19] proposed an active contour method based on feature learning, feature detection, and tracking using a Kalman filter. The semi-automatic tracking methods for the lumbar part of the spine are applied in research [16, 17] and in clinical practice [7, 20].

The methods published in literature for tracking the cervical vertebrae using DVF are limited. Reinartz et al. [15] proposed a method for tracking C0 to C6 based on template matching and normalized gradient field. Cervical vertebrae tracking is more complicated than tracking the lumbar vertebrae, because the former are smaller and more difficult to identify in the videos than the latter, due to the anatomy of the C1 and C2. C1 and C2 vertebrae have unique shapes, and thus they often appear as one single vertebra in DVF frames. What is more, the C7 vertebra is often partly shadowed by the shoulder. The aim of this study was to develop an efficient, semi-automatic method for tracking the cervical vertebrae in order to accurately calculate the intervertebral rotation. Our method required very little initial input and did not require an expertise in radiography. It was designed to continuously track the C1–C7 cervical vertebrae during flexion and extension movements, recorded in the sagittal plane using digital videofluoroscopy.

2 Method

The Phillips BV Libra mobile diagnostic X-ray image acquisition and viewing system was used to record the neck movements at a rate of 25 frames per second, with a frame size of 576 \(\times \) 720 pixels, and image resolution of 24 bits per pixel. The system consisted of two main components: (1) the C-arm containing the CCD camera, the image intensifier, the collimator and an X-ray tube; and (2) the mobile viewing station. Videos from two subjects were used for design and implementation of the method and videos from one subject was used for the proof of concept validation. For each subject two videos were recorded: one while the subject performed a flexion movement, and one during an extension movement.

The method we developed consisted of four steps: video preprocessing, manual input and template definition, tracking using template matching and Kalman filtering, and calculation of intervertebral rotation. We validated our method by comparing the results with a digitalized manual method using two fluoroscopic videos, where vertebral corners had been manually defined as landmarks by a trained clinician, which were used for calculation of intervertebral rotation based on DCRA.

2.1 Video Preprocessing

The fluoroscopic video frames were bandpass-filtered in order to enhance the contrast between the vertebral bodies and the surrounding tissue, and to remove the high frequency noise. The applied bandpass filter had a lower cut-off frequency of 5, and an upper cut-off frequency of 71, as proposed by Teyhen et al. [16]. Subsequently, a 10 \(\times \) 10 Canny edge detection filter with the lower threshold of 0.02 and the upper threshold of 0.05 was implemented to detect the edges of the vertebral bodies. Figure 1 shows each of the steps in the preprocessing.

Fig. 1.
figure 1

A segment of an original frame, a bandpass filtered frame and a fully preprocessed frame.

Full size image

2.2 Manual Input and Template Definition

A template for each of the cervical vertebrae was created by manual definition of landmarks in the first frame in the video, as illustrated in Fig. 2. These landmarks were used to define the templates in the first preprocessed frame and consequently each template was binary. Three separate methods were used for the template definition: one for C1, one for C2 and C7, and one for C3–C6.

Fig. 2.
figure 2

Illustration of the manually defined landmarks used for template definition in the first frame.

Full size image

2.2.1 C1

The template for C1 covered the vertebral arch which was visible in all anatomical positions during both flexion and extension movements. Landmarks were defined at the posterior and the anterior tube and at the posterior and anterior edges of the lateral masses to enable calculation of the intervertebral rotation, as illustrated in Fig. 2. The template was defined using the landmark at the posterior edge of the lateral masses and three additional landmarks positioned on the most cranial point, the most caudal point, and the most posterior point of the vertebral arch.

2.2.2 C2 and C7

The cranial part of the C2 vertebral body was not visible in all anatomical positions during movements due to shadowing by C1. C7 was also partially shadowed in the caudal part in some subjects by their shoulder. The corners of the vertebral bodies were defined for the computation of the intervertebral rotation and another four points were defined to enclose the template. The landmarks used for the definition of the template for C2 were positioned so the template enclosed the caudal part of the vertebral body and the cranial part of the vertebral arch. The anterior point was positioned just outside the anterior edge of the vertebral body and the cranial posterior point was positioned so that the template included the cranial part of the spinous process. For C7, the landmarks were positioned to enclose only the cranial part of the vertebral body.

2.2.3 C3–C6

The templates for C3 to C6 were defined using only the corners of the vertebral bodies. Four new points were automatically computed by adding eight pixels to each of the corners in the opposite direction of the diagonally placed corner. Those points were then connected, and this formed a quadrilateral, which constituted the template for the C3–C6 vertebrae. The size of the template was a compromise on having the template surrounding the entire vertebral body, but at the same time not containing other elements, such as a part of an adjacent vertebra.

2.3 Tracking Using Template Matching

The templates were tracked throughout the preprocessed video frames using cross-correlation and Kalman filtering. Cross-correlation was calculated for all combinations of rotation and spatial position of the templates in the frames within a limited search space. The search space was large enough to track the vertebrae movements, but small enough to avoid matching the template with an adjacent vertebra. The rotational search space was defined as the angle of the vertebra’s midplane in the previous frame ±4\(^{\circ }\), with a step size of 0.01\(^{\circ }\). The spatial search space was defined by adding 10 pixels to the position of the template in the previous frame in all directions. Subsequently, the maximum cross-correlation was found and used to determine the rotation and position resulting in the best match for the templates in the frame. Unrealistic high changes in the rotation and position of the templates between two adjacent frames were observed. Kalman filters have earlier been used to estimate the position of objects in images [18]. This inspired us to implement two Kalman filters in the tracking of the templates. One filter estimated the position of the templates by modelling the position and the velocity of the templates in both x and y dimensions. The other Kalman filter estimated the rotation by modelling the angle of the templates and their angular velocity. The acceleration was assumed to be constant for both the position and the rotation of the templates.

2.4 Intervertebral Rotation

The intervertebral rotation was calculated by means of the DCRA [11]. For the C2–C7 vertebrae, the corners of their vertebral bodies were used to calculate the anterior and posterior midpoints which defined midplanes of each vertebrae as illustrated in Fig. 3. The intervertebral rotation was calculated as the angle between the midplanes for two adjacent vertebrae. [11]. The corners of the C1 were detected in a different manner, as the C1 does not have a vertebral body. Instead, a midplane was defined as a line passing through the landmarks positioned at the anterior and the posterior tube (Fig. 2, upper left image).

Fig. 3.
figure 3

An illustration of the landmarks on a pair of adjacent vertebrae for the calculation of the intervertebral rotation. The corners of the vertebral bodies are marked with purple points. (Color figure online)

Full size image

2.5 Validation

We validated our semi-automatic method for vertebrae tracking through a comparison with a manual method using two fluoroscopic videos. The videos were recorded in the sagittal plane, and the subject performed movements of a full flexion and full extension. In the manual marking method, a trained clinician manually defined vertebral landmarks in 11 frames, constituting the 10% epochs of the total range of motion of movement. The midplanes for C3–C6 were defined in the same manner as in our semi-automatic method, whereas the midplanes for C1, C2, and C7 had been defined using only two landmarks. This entailed a modified version of the DCRA in the computation of the intervertebral rotation for these vertebrae in the manual method. The landmarks had been defined three times to enable examination and comparison of the variability both in our semi-automatic and in the manual method, which resulted with three datasets for each method. These three sets of manually defined landmarks in the first frame were used for the template definitions, thus the templates were defined three times. The intervertebral rotation was computed in all frames in the video and compared to the corresponding rotation computed by the manual method.

2.6 Intraobserver Repeatability

The Bland-Altman method was used to evaluate the intraobserver repeatability, and the differences between the intervertebral rotation were compared pairwise for both methods. The mean (d) and standard deviation (SD) were calculated for each set of differences between two datasets, and subsequently the Coefficient of Repeatability (CR) (Eq. 1), and the upper and lower Limits of Agreement (LOA) (Eq. 2) [4] were calculated as follows:

$$\begin{aligned} CR = 1.96 \cdot SD \end{aligned}$$
(1)
$$\begin{aligned} LOA = d \pm 1.96 \cdot SD \end{aligned}$$
(2)
Fig. 4.
figure 4

Intervertebral rotation for all pairs of adjacent cervical vertebrae during an extension (left column) and a flexion movement (right column). The red lines illustrate the intervertebral rotation obtained by our semi-automatic cervical spine tracking method for all frames, and the black asterisks illustrate the intervertebral rotation obtained with the manual method for the 11 manually marked frames. (Color figure online)

Full size image

3 Results

The tracking was qualitatively evaluated by means of visual inspection of the template tracking in the videos performed by our method. The evaluation did not reveal any tracking issues, except for minor irregularities in the end of the extension movement for two out of three C1 templates. The continuously intervertebral rotation both during extension and flexion, calculated by our semi-automatic method, is illustrated in Fig. 4 (red curves). The black asterisks illustrate the intervertebral rotation calculated from the manually marked vertebral landmarks by a trained clinician. The intervertebral rotation calculated by our method during an extension movement was comparable to the rotation calculated by the manual method for the C2/C3 to C6/C7 segments. However, differences were observed for the C1/C2, which might be due the irregularities in the tracking of the template for C1. For the flexion movement, the intervertebral rotation was comparable between the two methods for the C3/C4 to C6/C7 segments, but for C1/C2 and C2/C3 the rotation changed differently for the two methods throughout the movement. No issues were observed for the tracking of the templates for C2 and C3 and consequently, the observed differences might be a result of the midplane for C2 being differently defined in the manual method and our semi-automatic method.

Fig. 5.
figure 5

The Bland-Altman plot for the intervertebral rotation between the C1/C2 and C2/C3 segments.

Full size image

3.1 Intraobserver Analysis

The Bland-Altman plots were made for all pairs of adjacent vertebrae for intervertebral rotation for both our cervical vertebrae tracking method (abbreviated as CVTM in Fig. 5 and Tables 1 and 2), and the digital manual method, abbreviated as DMM. The Bland-Altman plots for intervertebral rotation between C1/C2 and C2/C3 are illustrated in Fig. 5. These plots illustrate the pair of adjacent vertebrae having the highest coefficient of repeatability. The Bland-Altman statistics for all pairs of adjacent vertebrae are presented in Table 1 for our semi-automatic approach and the manual method involving vertebrae marking by a trained clinician. In five out of six cases, the mean of differences for the CVTM was closer to zero than the mean of differences obtained from the DMM, indicating smaller bias in our method. The coefficient of repeatability for all pairs of adjacent vertebrae was lower for the semi-automatic method than for the manual method, except for the intervertebral rotation between C1/C2. This indicated that we obtained a higher intraobserver repeatability in most of the cases than if the intervertebral rotation was calculated using the manual marking approach.

Table 1. The Bland-Altman statistics for intervertebral rotation obtained by the CVTM and DMM for all pairs of adjacent vertebrae.
Full size table

4 Discussion

The purpose of this paper was to develop a semi-automatic tracking method for cervical vertebrae using template matching, in order to accurately and efficiently calculate intervertebral rotation of the cervical spine. The method enabled the tracking of C1 to C7 vertebrae in fluoroscopic videos during neck flexion and extension in the sagittal plane. Our semi-automatic method was validated as proof of concept through a comparison of the results with a digitalized manual method of vertebral marking performed by a trained clinician. The two methods showed comparable intervertebral kinematics. Additionally, our approach had a lower coefficient of repeatability than the manual method for all pairs of adjacent vertebrae, except for the intervertebral rotation of C1/C2. We acknowledge the need for a larger data set for validation, but the promising results obtained with our method indicate it to be a valid proof of concept for a fast and accurate calculation of intervertebral kinematics, with a potential of including intervertebral translation.

The user interaction with our semi-automatic is limited only to the initial step involving a definition of landmarks in the first video frame, making it far less time-consuming than the digitalized manual method. Our method had a low coefficient of repeatability, and accordingly was robust to small errors in the definition of landmarks for the templates, allowing successful template definition and intervertebral rotation calculations without an extensive expertise in radiography. The manual method of vertebrae marking requires a skilled clinician to define the landmarks in each frame to be analysed. Moreover, the manual method is susceptible to a variety of external factors, such as the consistency of clinician’s precision in marking the vertebrae.

No studies evaluating the repeatability for tracking of the cervical vertebrae were found in the literature. Reinartz et al. [15] tracked C0–C6 by a semi-automatic method and stated that the tracking performed by the algorithm was as accurate as manual template matching, however, no results for the validation of the method was presented. Thus, our method was compared to existing methods for calculation of intervertebral kinematics of the lumbar spine. Frobin et al. [11] calculated the intervertebral rotation between L3/L4 and L5/S1 in X-ray images. They defined the landmarks manually in all frames and obtained a CR of 1.94\(^\circ \). Teyhen et al. [16] calculated intervertebral rotation for L3/L4 to L5/S1 using a semi-automatic method for DVF. The method was developed to optimize the precision of manually defined landmarks in all frames, and it resulted in a CR of 2.31\(^\circ \). Penning et al. [14] proposed a method based on image registration for calculation of intervertebral rotation between L1/L2 to L5/S1 in X-ray images and this method had a CR at 0.86–0.94\(^\circ \). Yeager et al. [20] evaluated a vertebral tracking algorithm, called KineGraph Vertebral Motion Analysis (VMA), which is used to calculate intervertebral rotation between L1/L2 to L5/S1 in DVF and the coefficient of repeatability of this method was 1.53\(^\circ \). An overview of the repeatability coefficients is presented in Table 2, where we compared our method to the aforementioned studies found in literature. However, methods for tracking lumbar and cervical vertebrae cannot be directly compared, since the cervical vertebrae are smaller than the lumbar vertebrae. Consequently, the intervertebral rotation for the cervical vertebrae are more sensitive to variability in the definition of the landmarks and thereby it is more difficult to obtain a low CR.

Table 2. CR and methods for calculation of intervertebral rotation for the CVTM and the methods found in the literature.
Full size table

The semi-automatic method for cervical vertebrae tracking presented in our study had a lower CR for the C2/C3 to C6/C7 segments than Frobin et al. [11], Teyhen et al. [16], and Yeager et al. [20]. The low CR for our semi-automatic method compared to other methods found in literature might be due to our implementation of the Kalman filters in vertebrae tracking. The Kalman filters estimate a more precise position and rotation of the templates. The high CR for the intervertebral C1/C2 rotation is primarily a result of irregularities in the tracking of the template for C1 during the extension movement. The irregularity in the C1 tracking was likely caused by a change in the shape of C1 during the movement, which resulted in a deviation from the assumption of rigidity. This may be due to a rotation around the longitudinal axis. Consequently, a differently defined template may improve the tracking of C1 and thereby decrease the CR for the intervertebral rotation between C1/C2.