Robust vertebra identification using simultaneous node and edge
  predicting Graph Neural Networks
Automatic vertebra localization and identification in CT scans is important
for numerous clinical applications. Much progress has been made on this topic,
but it mostly targets positional localization of vertebrae, ignoring their
orientation. Additionally, most methods employ heuristics in their pipeline
that can be sensitive in real clinical images which tend to contain
abnormalities. We introduce a simple pipeline that employs a standard
prediction with a U-Net, followed by a single graph neural network to associate
and classify vertebrae with full orientation. To test our method, we introduce
a new vertebra dataset that also contains pedicle detections that are
associated with vertebra bodies, creating a more challenging landmark
prediction, association and classification task. Our method is able to
accurately associate the correct body and pedicle landmarks, ignore false
positives and classify vertebrae in a simple, fully trainable pipeline avoiding
application-specific heuristics. We show our method outperforms traditional
approaches such as Hungarian Matching and Hidden Markov Models. We also show
competitive performance on the standard VerSe challenge body identification
task.
Vincent Bürgin, Raphael Prevost, Marijn F. Stollenga
# Robust verteb identification using simultaneous node and edge predicting Graph Neural Networks

###### Abstract

Automatic vertebra localization and identification in CT scans is important for numerous clinical applications. Much progress has been made on this topic, but it mostly targets positional localization of vertebrae, ignoring their orientation. Additionally, most methods employ heuristics in their pipeline that can be sensitive in real clinical images which tend to contain abnormalities. We introduce a simple pipeline that employs a standard prediction with a U-Net, followed by a single graph neural network to associate and classify vertebrae with full orientation. To test our method, we introduce a new vertebra dataset that also contains pedicle detections that are associated with verteb bodies, creating a more challenging landmark prediction, association and classification task. Our method is able to accurately associate the correct body and pedicle landmarks, ignore false positives and classify vertebrae in a simple, fully trainable pipeline avoiding application-specific heuristics. We show our method outperforms traditional approaches such as Hungarian Matching and Hidden Markov Models. We also show competitive performance on the standard VerSe challenge body identification task.

Keywords:graph neural networks spine localization spine classification deep learning.

## 1 Introduction

Vertebra localization and identification from CT scans is an essential step in medical applications, such as pathology diagnosis, surgical planning, and outcome assessment [6, 3]. This is however a tedious manual task in a time-sensitive setting that can benefit a lot from automation. However, automatically identifying and determining the location and orientation of each vertebra from a CT scan can be very challenging: (i) scans vary greatly in intensity and constrast, (ii) metal implants and other materials can affect the scan quality, (iii) vertebrae might be deformed, crushed or merged together due to medical conditions, (iv) vertebrae might be missing due to accidents or previous surgical operations.

Recently, public challenges like the VerSe challenge [24] have offered a common benchmarking platform to evaluate algorithms to automate this task, resulting in a boost in research on this topic. However, these challenges focus on
finding the position of vertebrae, ignoring the orientation or direction. Additionally, practically all methods employ manual heuristic methods to identify landmarks and filter out false positives.

In this paper, we introduce a trainable method that performs vertebrae localization, orientation estimation and classification with a single architecture. We replace all hand-crafted rules and post-processing steps with a single trainable Graph Neural Network (GNN) that learns to filter out, associate and classify landmarks. We apply a generalized Message Passing layer that can perform edge and node classification simultaneously. This alleviates the need for sensitive hand-tuned parameter tuning, and increases robustness of the overall pipeline.

The main contributions of our work are: (1) introducing a pipeline that uses a single Graph Neural Network to perform simultaneous vertebra identification, landmark association, and false positive pruning, without the need for any heuristic methods and (2) building and releasing a new spine detection dataset that adds pedicles of vertebrae to create a more complex task that includes orientation estimation of vertebrae, which is relevant for clinical applications.

## 2 Related Work

**Spine Landmark Prediction and Classification** The introduction of standardised spine localization and classification challenges [10, 24] resulted in a boost in research on this problem. Convolutional Neural Networks (CNN) became a dominant step in most approaches shortly after their introduction [7]. Most modern methods process the results of a CNN using a heuristic method to create a 1-dimensional sequence of vertebra detections, before applying classification: [27] generate heatmaps for body positions, and refine it into a single sequence graph that uses message passing for classification. [14] generate heatmaps, extract a 1-dimensional sequence and use a recurrent neural network for classification. [18] produce a heatmap using a U-Net [21], but use a simpler approach by taking the local maxima as landmarks, and forming a sequence by accepting the closest

Figure 1: Summary of the proposed method. A CNN generates a heatmap from which local maxima are connected using k-nearest neighbours to form a graph. Then a single GNN associates the keypoints, performs classification and filters out false postives. The full pipeline is trained and does not require hand-tuned post-processing.
vertebra that is within a defined distance range. [16] uses a directional graph and Dynamic Programming to find an optimal classification.

**Graph Neural Networks** In recent years, Graph Neural Networks (GNN) have surged in popularity [11], with a wide and growing range of applications [26]. A prominent task in the literature is node-level representation learning and classification. A less prominent task is edge classification, for which early work used a dual representation to turn edge- into node representations [1]. Other approaches model edge embeddings explicitly, such as [12] and [4]. The most general formulation of GNNs is the message-passing formulation [5] which can be adapted to perform both edge and node classification at the same time. We use this formulation in our method.

Various methods have applied GNNs to keypoint detection, however they all apply to 2-dimensional input data. In [20] GNNs are used to detect cars in images. An edge classification task is used to predict occluded parts of the car. However, the GNN step is ran individually for every car detection and the relation between cars is not taken into account, unlike our task. Also there is no node classification applied. In [15] a GNN is used to group detected keypoints for human-pose estimation on images. Keypoints are grouped using edge prediction where edge-embeddings are used as input the the GNN. A separate GNN processes node embeddings to facilitate the final grouping. However, the node and edge embeddings are processed separately from each other.

## 3 Method

In this paper we tackle vertebra localization and classification, but unlike other methods that only focus on detecting the body of the vertebrae, we also detect the pedicles and associate them with their corresponding body. This allows us to also define the _orientation_ of the vertebra defined by the plane passing through the body and pedicles 1. To this end we introduce a new dataset that includes pedicles for each vertebra, described in Section 4.1, creating a challenging keypoint detection and association task.

Footnote 1: Another common way to define the orientation is using the end-plates of the vertebra body; however end-plates can be irregular and ill-defined in pathological cases

Our method consists of a two-stage pipeline shown in Figure 1: we detect keypoints from image data using a _CNN stage_, and form a connected graph from these keypoints that is processed by a _GNN stage_ to perform simultaneous node and edge classification, tackling classification, body to pedicle association as well as false positive detection with a single trainable pipeline without heuristics.

**CNN stage** The CNN stage detects candidate body, left pedicle and right pedicle keypoints and provides segment classifications for the body keypoints as either cervical, thoracic, lumbar or sacral. We use a UNet2 CNN [19] and select all local maxima with an intensity above a certain threshold \(\tau\), in this paper we use \(\tau=0.5\). These keypoints are connected to their \(k\) nearest neighbours, forming a graph. In rare cases this can result in unconnected cliques, in which case the nearest keypoint of each clique is connected to \(\frac{k}{3}\) nearest points in
the other cliques, ensuring a fully connected graph. All nodes and edges are associated with information through embeddings, described below.

**GNN stage** The second stage employs a generalized message-passing GNN following [5] to perform several prediction tasks on this graph simultaneously:

1. **keypoint association prediction:** we model association between body keypoints and their corresponding pedicle keypoints as binary edge classification on the over-connected \(k\)-NN graph.
2. **body keypoint level prediction:** for body keypoints, we model the spine level prediction as multi-class node classification.
3. **keypoint legitimacy prediction:** to filter out false-positive keypoints, we additionally compute an binary legitimacy prediction for each node.

To perform these task, our message-passing GNN maintains edge and node embeddings which are updated in each layer. A message-passing layer performs a node update and edge update operation. Denoting the feature vector of a node \(v\) by \(x_{v}\), and the feature vector of a directed edge \((u,v)\) by \(x_{uv}\), the node and edge features are updated as follows:

\[\underbrace{x_{u}^{\prime}=\bigoplus_{v\in\mathcal{N}_{u}\cup\{u\}}\psi_{ \text{node}}(x_{u},x_{v},x_{uv})}_{\text{Node update}},\qquad\underbrace{x_{uv}^{ \prime}=\psi_{\text{edge}}(x_{u},x_{v},x_{uv})}_{\text{Edge update}} \tag{1}\]

Here \(\bigoplus\) denotes a symmetric pooling operation (in our case max pooling) over the neighborhood \(\mathcal{N}_{u}\). \(\psi_{\text{node}}\) and \({}_{\text{edge}}\) are trainable parametric functions: in our case, two distinct two-layer MLPs with ReLU nonlinearities. After \(N\) such message-passing layers we obtain an embedding vector for each node and edge. Each node/edge embedding is passed through a linear layer (distinct for nodes and edges) to obtain a vector of node class logits or a single edge prediction logit, respectively. The last entry in the node prediction vector is interpreted as a node legitimacy prediction score: nodes predicted as illegitimate are discarded for the output.

The node input features \(x_{u}\in\mathbb{R}^{7}\) consist of the one-hot encoded keypoint type (body, left or right pedicle) and the segment input information (a pseudo-probability in \([0,1]\) for each of the four spine segments of belonging to that segment, computed by applying a sigmoid to the heatmap network's output channels which represent the different spine segments). The edge input features \(x_{uv}\in\mathbb{R}^{4}\) consist of the normalized direction vector of the edge and the distance between the two endpoints.

The output of the GNN contains finer spine-level classification (i.e. C1-C7, T1-T13, L1-L6, S1-S2), keypoint-level legitimacy (legitimate vs. false-positive detection) and body-pedicle association via edge prediction, implicitly defining the orientation of each vertebra. Prediction scores of corresponding directed edges \((u,v)\) and \((v,u)\) are symmetrized by taking the mean.

In our experiments we consider variations to our architecture: weight sharing between consecutive GNN layers, multiple heads with a shared backbone (jointly trained) and dedicated networks (separately trained) for edge/node prediction.
## 4 Experiments

### Dataset

Our main dataset consists of 2118 scans, of which 1949 form the training dataset, and 169 the validation dataset. This includes 360 scans from the VerSe datasets [24], 1676 scans from the CT Colonography dataset [25] and 82 scans from the CT Pancreas dataset [22]. Of these datasets, only VerSe has labels for spine levels, and none of the datasets have labels for the pedicles of vertebrae. The labeling was done in a two-step process with initial pedicle locations determined through post-processing of ground truth segmentations, then transfered to other datasets via bootstrapping. For each verteb, the verteb level is labeled, including the position of the body and the right and left pedicles. This keypoint dataset is publicly available at [https://github.com/ImFusionGmbH/VID-vertebra-identification-dataset](https://github.com/ImFusionGmbH/VID-vertebra-identification-dataset).

Additionally, we also perform the VerSe body classification task on the original VerSe dataset [24] which contains 160 scans in an 80/40/40 split.

### Training

Heatmap networkThe heatmap is generated by a UNet\({}^{2}\) network [19], with 4 layers and 32 channels. The network is trained on crops of \(128^{3}\) voxels with 1.5mm spacing. The target data consists of Gaussian blobs (\(\sigma=6\)mm) at the landmark positions. We sample 70% of the crops around landmark locations, and 30% randomly from the volume. Additionally we apply 50% mirror augmentations, and 20 degree rotations around all axes. We use the MADGRAD optimizer [8] with a learning rate of \(10^{-3}\). We use a binary cross-entropy loss, with an 80% weighting towards positive outputs to counter the data's balancing.

Graph Neural NetworkThe GNN is implemented in _PyTorch Geometric_2.0.4 [9]. The three predictions of the graph neural network - edge classification, node classification and node legitimacy prediction - are trained via cross-entropy losses which are weighted to obtain the overall loss:

\[\mathcal{L}=\alpha\mathcal{L}_{\text{edge}}^{\text{BCE}}+\beta\mathcal{L}_{ \text{node}_{\text{class}}}^{\text{CE}}+\gamma\mathcal{L}_{\text{node}_{ \text{legit}}}^{\text{BCE}} \tag{2}\]

Figure 2: Examples of challenging cases for the Hungarian Matching edge detection baseline, and corresponding correct detections by our GNN architecture (red = false positive, blue = false negative detection).
We only make edge predictions on edges that run between a body and a pedicle keypoint - the other edges are only used as propagation edges for the GNN. Similarly, we only make spine level predictions on body keypoints. Solely these subsets of nodes/edges go into the respective losses.

As input data, we use the predicted keypoints of the heatmap network on the training/validation set. The ground-truth keypoints are associated to this graph to create the target of the GNN. We include three synthetic model spines during training (keypoints in a line spaced 30mm apart) to show the network typical configurations and all potential levels (with all levels/without T13, L6, S2/without T12, T13, L6, S2).

We tune various hyperparameters of our method, such as network depth and weight sharing, the \(k\) of the \(k\)-NN graph, the loss weighting parameters \(\alpha,\beta,\gamma\) and the number of hidden channels in the message-passing MLP. We use a short notation for our architectures, such as (5x1, 4, 1) for 5 independent message-passing layers followed by 4 message-passing layers with shared weights and another independent message-passing layer.

Various data augmentations are used to make our network more robust to overfitting: (i) rotation of the spine by small random angles, (ii) mirroring along the saggital axis (relabeling left/right pedicles to keep consistency), (iii) perturbation of keypoints by small random distances, (iv) keypoint duplication and displacement by a small distance (to emulate false-positive duplicate detections of the same keypoint), (v) keypoint duplication and displacement by a large distance (to emulate false-positive detections in unrelated parts of the scan) and (vi) random falsification of spine segment input features. We define four levels of augmentation strength (no/light/default/heavy augmentations) and refer to the supplementary material for precise definitions of these levels.

### Evaluation

We evaluate our method on two tasks. The first one is the full task consisting of node and edge classification for vertebra level and keypoint association detection. We evaluate this on the 2118 scan dataset which comes with pedicle annotations. The second task consists only of vertebra localization, which we evaluate on the VerSe 2019 dataset [24]. Unless otherwise specified, we use \(k=14\), the (13x1) architecture, batch size 25 and reapparent every 25 epochs for the full task, and \(k=4\), the (9x1) architecture, batch size 1 and reapparent every epoch for the VerSe task. In both cases we use the default augmentation level and \(\alpha=\beta=1\). As evaluation metrics we use the VerSe metrics _identification rate_ (ratio of ground-truth body keypoints for which the closest predicted point is correctly classified and within 20mm) and \(d_{mean}\) (mean distance of correctly identified body keypoints to their 1-NN predictions) [24]. Furthermore we evaluate the edge and illegitimacy binary predictions by their \(F_{1}\) scores. Since the identification rate is largely unaffected by false-positive predicted keypoints, we disable legitimacy predictions unless for specific legitimacy prediction experiments to help comparability.
We compare our methods to two baselines: For node prediction, we use a Hidden Markov Model (HMM) [2] that is fitted to the training data using the Baum-Welch algorithm using the _pomegranate_ library [23]. The HMM gets the predicted segment labels in sequence. Dealing with false-positive detection outliers is very difficult for this baseline, therefore we filter out non-legitimate detections for the HMM inputs to get a fairer comparison, making the task slightly easier for the HMM. For the VerSe challenge, we also compare our results to the top papers reported in the VerSe challenge [24]. For edge prediction, we compare our method to Hungarian matching on the keypoints from the CNN pipeline.

### Results

The results on our introduced dataset are shown in Table 1. Our method outperforms the baseline methods, both in identification rate and edge accuracy. We also show results on a 'hard' subset, selected as the data samples where either of the methods disagree with the ground truth. This subset contains 47 scans and represents harder cases where the association is not trivial. We give \(p\)-values of the Wilcoxon signed-ranked test against the respective baseline for numbers that are better than the baseline. Figure 2 shows a qualitative example where the baseline method fails due to false-positive and false-negative misdetections in the input (these errors in the input were generated by augmentations and are more extreme than usual, to demonstrate several typical baseline failures in one example). The GNN correctly learns to find the correct association and is not derailed by the misdetections.

The examples where our architecture fails typically have off-by-one errors in the output from the CNN: for example, the last thoracic vertebra is detected as a lumbar vertebra (usually in edge cases where the two types are hard to distinguish). Hence the GNN classifications of the lumbar segment, and possibly the thoracic segment, will be off by one (see Figure S2 in the supplementary).

Table 2 shows the performance differences of various GNN architectures. The single-head architecture with 13 individual layers performs the best on identi

\begin{table}
\begin{tabular}{c c|c c c} \hline \hline  & Method & identification rate & edge \(F_{1}\) score & \(d_{\text{mean}}\) \\ \hline \multirow{4}{*}{Edge vs. node classification} & Joint node/edge & \(\begin{array}{c}\textbf{97.19}\,/\,\textbf{89.88}\\ (p\,=\,0.080)\end{array}\) & \(\begin{array}{c}\textbf{98.81}\,/\,\textbf{96.22}\\ \end{array}\) & \(\textbf{1.68}\,/\,\textbf{1.95}\) \\  & Node only & \(\begin{array}{c}\textbf{96.91}\,/\,\textbf{88.90}\\ (p\,=\,0.109)\end{array}\) & — & \(\begin{array}{c}\textbf{1.68}\,/\,\textbf{1.96}\\ \end{array}\) \\  & Edge only & — & \(\begin{array}{c}\textbf{99.31}\,/\,\textbf{97.77}\\ (p\,=\,0.019)\end{array}\) & — \\ \hline \multirow{2}{*}{Baselines} & Hidden Markov & \(\begin{array}{c}\textbf{94.29}\,/\,\textbf{79.46}\\ \end{array}\) & — & \(\begin{array}{c}\textbf{1.81}\,/\,\textbf{2.52}\\ \end{array}\) \\  & Hungarian matching & — & \(\begin{array}{c}\textbf{98.93}\,/\,\textbf{96.56}\\ \end{array}\) & — \\ \hline \hline \end{tabular}
\end{table}
Table 1: Results on the 2118 spine dataset (full validation set / hard subset). Comparing best GNN architectures with the baselines. Wilcoxon signed-rank test \(p\)-values are given for GNN numbers that outperform the baseline (by the test’s construction, \(p\)-values agree between the full and hard subset).
fication rate, although the 13-layer architecture with 11 shared layers performs very similarly. The architectures with fewer layers perform slightly better in edge accuracy since this task is less context dependent. Multi-head architectures perform slightly worse on identification, but retain a good edge accuracy. Training a model solely directed to edge classification does yield the highest edge accuracy, as expected. Enabling legitimacy predictions slightly degrades the performance, likely due to two facts: for one, an extra loss term is added, which makes the

\begin{table}
\begin{tabular}{c c|c c c} \hline \hline \multicolumn{1}{c|}{Method} & identification rate & edge \(F_{1}\) score & illegitimacy \(F_{1}\) \\ \hline \multirow{4}{*}{Single-head architectures} & (7x1) & 96.74 & 98.86 & — \\  & (13x1) & **97.19** & 98.81 & — \\  & (1,5,1) & 96.75 & 98.91 & — \\  & (1,11,1) & 97.10 & 98.67 & — \\ \hline Multi-head architectures: backbone, edge/node head & (1x1), (4x1), (12x1), (10x1) & 96.94 & 99.16 & — \\  & (3x1), (2x1), (10x1) & 96.87 & 99.00 & — \\  & (5x1), (—), (8x1) & 96.79 & 98.88 & — \\ \hline \multirow{4}{*}{Dedicated edge architectures} & (3x1) & — & 99.28 & — \\  & (5x1) & — & 99.28 & — \\  & (1,3,1) & — & **99.35** & — \\ \hline \multirow{4}{*}{With legitimacy prediction: different logit. loss weights} & \(\gamma=0.1\) & 95.61 & **98.96** & 56.00 \\  & \(\gamma=1.0\) & **96.50** & 98.75 & 62.41 \\  & \(\gamma=5.0\) & 96.35 & 98.18 & 61.67 \\  & \(\gamma=10.0\) & 96.42 & 98.96 & **62.78** \\ \hline \multirow{4}{*}{Augmentations} & None & 95.27 & 96.83 & — \\  & Light & 97.18 & **98.83** & — \\  & Default & **97.19** & 98.81 & — \\  & Heavy & 97.03 & 98.67 & — \\ \hline \hline \end{tabular}
\end{table}
Table 2: Hyperparameter comparisons on the 2118 spine dataset. Comparing GNN architectures, influence of enabling legitimacy predictions and of different augmentation strengths. Best values within each group highlighted in bold.

\begin{table}
\begin{tabular}{c c|c c} \hline \hline \multicolumn{1}{c|}{Method} & identification rate & \(d_{\text{mean}}\) & illegitimacy \(F_{1}\) \\ \hline Single-head GNN & (9x1) & \(\begin{array}{c}93.26\,/\,93.02\\ (p=3.1e-6/p=5.2e-7)\end{array}\) & \(\begin{array}{c}1.28\,/\,1.43\\ \end{array}\) & — \\ \hline With legitimacy prediction & \(\gamma=10.0\) & \(\begin{array}{c}90.75\,/\,87.51\\ (p=4.8e-6/p=9.4e-7)\end{array}\) & \(\begin{array}{c}\textbf{1.23}\,/\,\textbf{1.32}\\ \end{array}\) & 77.67\(\,/\,81.69\) \\ \hline Baseline & Hidden Markov & \(48.59\,/\,49.06\) & \(1.32\,/\,1.45\) & — \\ \hline \multirow{2}{*}{VerSe challenge entries [24]} & Payer C. [17] & \(95.65\,/\,\textbf{94.25}\) & \(4.27\,/\,4.80\) & — \\  & Lessmann N. [13] & \(89.86\,/\,90.42\) & \(14.12\,/\,7.04\) & — \\  & Chen M. [24] & **96.94**\(\,/\,86.73\) & \(4.43\,/\,7.13\) & — \\ \hline \hline \end{tabular}
\end{table}
Table 3: Results on the VerSe 2019 dataset (validation/test set). We compare our GNN architecture, our Hidden Markov baseline, and the reported numbers of the three top VerSe challenge entries [24]. Best values highlighted in bold.
model harder to train. Also, the identification rate metric is not majorly affected by having additional false-positive detections, hence there is little to be gained in terms of this metric by filtering out false positives. An optimal legitimacy loss weighting seems to be \(\lambda=1.0\). Finally, augmentations help generalization of the model, but the amount of augmentations seems to have little effect.

Table 3 shows the result on the traditional VerSe body-identification task. Our method yields a competitive performance despite not being optimized on this task (identification rate slightly lower than the leaders but a better average distance to the landmarks).

## 5 Conclusion

We introduced a simple pipeline consisting of a CNN followed by a single GNN to perform complex verteba localization, identification and keypoint association. We introduced a new more complex verteba detection dataset that includes associated pedicles defining the full orientation of each verteba, to test our method. We show that our method can learn to associate and classify correctly with a single GNN that performs simultaneous edge and node classification.

The method is fully trainable and avoids most heuristics of other methods. We also show competitive performance on the VerSe body-identification dataset, a dataset the method was not optimized for. We believe this method is general enough to be usable for many other detection and association tasks, which we will explore in the future.
