Differently from state-of-the-art methods that use heatmap based regression for predicting single landmark locations, in this work we show that heatmaps can also be used to estimate the landmark uncertainty caused by annotation ambiguities. To this end, we first adapt the heatmap regression framework to allow optimization of not only the model parameters of the network but also the shape parameters of the target heatmaps. By doing so, we model the uncertainty based on annotation ambiguities of the entire training dataset, representing a homoscedastic aleatoric uncertainty that is the same for each image in the dataset. We name this uncertainty dataset-based annotation uncertainty. Second, during inference, a Gaussian function is fitted to the predicted heatmap to not only predict the most probable landmark location but its distribution, i.e., a heteroscedastic aleatoric uncertainty that is specific for each image. We name this uncertainty sample-based annotation uncertainty. We evaluate our approach on two datasets, i.e., on left hand radiographs and lateral cephalograms, showing that the uncertainties correlate with the magnitude and direction of the average localization error as well as the inter-observer variabilities.

This work extends our MICCAI-UNSURE 2020 workshop paper (Payer et al., 2020) in the following: (1) We improved the representativeness of the inter-observer variability used in our experiments by extending the annotation of the cephalogram dataset with additional annotations from nine observers. (2) We performed an in-depth analysis of our proposed method for estimating the dataset-based and sample-based annotation uncertainty on the newly annotated dataset. (3) We introduced a new experiment showing the importance of utilizing landmark location uncertainty into decision making with an example of measuring anatomical abnormalities in lateral cephalograms. (4) Finally, we included more details of the method, and extended our literature review and discussion section.

Due to the variable visual appearance of the anatomy defining a landmark position in images, modern methods for landmark localization are based on convolutional neural networks (CNNs) because of their superior image feature extraction capabilities as compared to other machine learning techniques (Toshev and Szegedy, 2014; Tompson et al., 2014; Wei et al., 2016; Payer et al., 2016; Newell et al., 2016; Yang et al., 2017; Payer et al., 2019; Bier et al., 2019). While earlier work that applied CNNs for human pose estmation – a task similar to anatomical landmark localization – uses CNNs for directly regressing landmark coordinates (Toshev and Szegedy, 2014), more recent methods for human pose estimation apply CNNs with the heatmap regression framework (Tompson et al., 2014; Wei et al., 2016; Newell et al., 2016). In the heatmap regression framework, the network is trained not to directly predict the landmark coordinates, but a heatmap which is generated as an image of a Gaussian function centered at the coordinate of the landmark. Thus, in this image-to-image mapping, the network learns to generate high responses on locations close to the target landmarks, while responses on wrong locations are being suppressed.

In our previous work (Payer et al., 2016), we introduced the CNN-based heatmap regression framework to anatomical landmark localization, while additionally proposing a CNN architecture called SpatialConfiguration-Net (SCN), which integrates spatial information of landmarks directly into an end-to-end trained, fully convolutional network. Building upon our approach, Yang et al. (2017) used the heatmap regression framework to generate predictions for landmarks with missing responses and incorporated an additional postprocessing step to remove false-positive responses. Bier et al. (2018) used heatmap regression for localizing landmarks on radiographs of the hips. To cope with the small number of images of their training dataset, they implement a sophisticated image synthesis framework based on projecting three-dimensional volumes to two-dimensional images and pretrain their networks with these synthetic images. Mader et al. (2019) adapted the heatmap regression networks for predicting more than a hundred landmarks on a dataset of spinal CT volumes. However, due to memory restrictions resulting from a large number of landmarks, they do not incorporate upsampling layers, thus reducing the accuracy of their heatmap regression outputs. Other recent methods use a coarse-to-fine approach by combining global and local image information with an attention-guided mechanism (Zhong et al., 2019), or with multiplying intermediate network outputs to restrict local responses to their global position (Chen et al., 2019).

In all aforementioned methods, the size of the target heatmaps, i.e., the variance of the Gaussian function, is fixed for each landmark and needs to be predefined before network training. Although our previous work in Payer et al. (2019) indicated that learning the parameters of an isotropic target Gaussian function for each landmark heatmap independently increases prediction accuracy, we did not provide any interpretation of the estimated heatmap size. Moreover, while there exist work that fits isotropic Gaussian functions to the predicted heatmap outputs for obtaining a more robust point prediction, e.g., Liu et al. (2019), most methods only take the single point-wise coordinate prediction as a final output and are not considering the remaining information in the heatmap prediction. Thus, to the best of our knowledge, none of the methods has investigated how the target Gaussian functions or the Gaussian function fitted to the heatmap prediction correlates with the uncertainty of landmark localization.

Uncertainty estimation in medical imaging has been widely adopted and used in segmentation, e.g., for segmenting multiple sclerosis lesions (Nair et al., 2020), the whole-brain (Roy et al., 2018), the fetal brain and brain tumors (Wang et al., 2019a), the carotid artery (Camarasa et al., 2020) as well as for ischemic stroke lesion and vessel extraction from retinal images (Kwon et al., 2020). Furthermore, uncertainty estimation was applied for disease detection using a diabetic retinopathy dataset (Leibig et al., 2017; Ayhan and Berens, 2018), for classification of skin lesions (Mobiny et al., 2019) and for probabilistic classifier-based image registration (Sedghi et al., 2019). A number of uncertainty metrics have been proposed in literature, which can be categorized into model (epistemic) and data (aleatoric) uncertainty (Kendall and Gal, 2017). While epistemic uncertainty represents the prediction confidence of the model for an unseen sample based on its similarity to training data samples, aleatoric uncertainty is caused by the observation noise inherent to the dataset.

Epistemic uncertainty for neural networks involves the use of probabilistic models, where the prediction given an image is not deterministic. This can be achieved with Bayesian Neural Networks (BNN), where weights are not represented as scalars but as distributions. During training, the parameters of the weight distributions are learned, while during inference, the probabilistic output is generated by combining multiple predictions for different weights sampled from the learned distributions (Blundell et al., 2015; Shridhar et al., 2019). A popular approximation of BNNs is Monte Carlo-Dropout (MCD) (Gal and Ghahramani, 2016). Here, the weight distributions are approximated by reusing Dropout (Srivastava et al., 2014) during inference to retrieve a set of stochastic predictions. Since MCD is easy to implement and computationally more efficient than BNNs, it is the most common method to model epistemic uncertainty in medical imaging (Leibig et al., 2017; Roy et al., 2018; Mobiny et al., 2019; Wang et al., 2019a; Camarasa et al., 2020; Kwon et al., 2020; Nair et al., 2020).

Aleatoric uncertainty for neural networks can be estimated by considering network outputs as distributions rather than point predictions. For classification tasks, this can be achieved by using not only the most probable class as the prediction, but the distribution represented by the probabilistic network outputs (Kendall and Gal, 2017). For regression tasks, this can be achieved by adding additional network outputs that model the assumed distribution of the observed annotation noise (Kendall and Gal, 2017; Kwon et al., 2020). For example, when assuming the annotation noise is normally distributed, the network is trained to not only predict the mean but also the variance of the distribution as an additional network output. Aleatoric uncertainty has only recently been investigated in landmark localization for computer vision tasks like pose estimation (Chen et al., 2019; Gundavarapu et al., 2019). Nevertheless, it was only done in the context of occluded landmarks. Due to large inter- and intra-observer variabilities (Chotzoglou and Kainz, 2019; Jungo et al., 2018), the dominant source of aleatoric uncertainty in the medical domain is annotation ambiguity. To the best of our knowledge, there is no previous work that shows how to evaluate annotation uncertainty in medical images from noisy data annotation.

In landmark localization the task is to use an input image ${I:\Omega_{I}\rightarrow\mathbb{R}}$ with $\Omega_{I}\subset\mathbb{R}^{2}$ to determine the coordinates $\mathbf{x}_{i}\in\Omega_{I}$ of $N$ specific landmarks $L_{i}$ within the extent $\Omega_{I}$ of the image $I$. In the heatmap regression framework, a CNN with parameters $\mathbf{w}$ is set up to predict $N$ heatmaps $\hat{h}_{i}:\Omega_{I}\rightarrow\mathbb{R}$, representing the pseudo-probability of the location of landmark $L_{i}$. As has been done in several other works (e.g., Tompson et al. (2014); Newell et al. (2016); Payer et al. (2016); Zhong et al. (2019)), a target heatmap ${h_{i}:\Omega_{I}\rightarrow\mathbb{R}}$ for landmark $L_{i}$ with coordinate $\mathbf{x}_{i}$ can be represented as an isotropic two-dimensional Gaussian function, i.e.,

The mean of the Gaussian is set to the target landmark’s coordinate ${\mathbf{x}_{i}\in\Omega_{I}}$, while the size of the isotropic Gaussian is defined by $\sigma_{i}$. A scaling factor $\gamma$ is used to avoid numerical instabilities during training. Fully convolutional networks that predict $N$ heatmaps $\hat{h}_{i}$ can then be trained with a pixel-wise $L_{2}$ loss to minimize the difference between predicted heatmaps $\hat{h}_{i}$ and target heatmaps $h_{i}$ for all landmarks $L_{i}$ by updating the network parameters $\mathbf{w}$, i.e.,

Crucial hyperparameters when optimizing the CNN are the sizes $\sigma_{i}$ of the Gaussian heatmaps, since suboptimal $\sigma_{i}$ lead to either reduced prediction accuracy or reduced robustness (Payer et al., 2019). All landmark specific $\sigma_{i}$ can be found via hyperparameter searching techniques (e.g., grid search), however, such approaches are time consuming. Therefore, in our previous work (Payer et al., 2019), sigmas which are unique for each landmark are learned directly during network optimization. Thus, in Payer et al. (2019), we adapted the pixel-wise loss function from Eq. (2) to

Since the $\sigma_{i}$ parameters of the Gaussian functions are also treated as unknowns, we enable learning them in addition to the network parameters $\mathbf{w}$. Thus, the gradient is not only backpropagated through the predicted heatmaps $\hat{h}_{i}$ to update the network parameters $\mathbf{w}$, but also through the target heatmaps $h_{i}$ to update the individual target heatmap sizes $\sigma_{i}$. The second term in Eq. (3) is used as a regularization to prefer smaller sigma values, as otherwise minimizing the difference between $\hat{h}_{i}$ and $h_{i}$ could lead to the trivial solution where $\sigma_{i}\rightarrow\infty$ with $h_{i}\approx 0$. It is important to note that although the learned target heatmap sizes are different for each landmark, each target heatmap size $\sigma_{i}$ is the same for all images in the dataset.

As we will show in our experiments, learning $\sigma_{i}$ values enables the network to encode the prediction uncertainty for each landmark individually. The more certain the network is in its prediction of a specific landmark, the smaller the target heatmap size will be. However, when using isotropic Gaussian functions as in Eq. (1) as the target heatmaps, information of the directional uncertainty cannot be modeled, which would be beneficial for landmarks that are ambiguous only in a specific direction (see Fig. 1).

In this work, we propose to model the target heatmaps with anisotropic Gaussian functions. A target heatmap is represented as an image of an anisotropic two-dimensional Gaussian function ${h_{i}(\mathbf{x};\mathbf{\Sigma}_{i}):\Omega_{I}\rightarrow\mathbb{R}}$, i.e.,

Differently to Eq. (1), the shape of the Gaussian is not defined by only a single $\sigma_{i}$ value but by a full covariance matrix $\mathbf{\Sigma}_{i}$. This way, the target heatmaps are able to have an orientation and varying extents in different directions. To separately analyze the orientation and the extent, the covariance matrix $\mathbf{\Sigma}_{i}$ can be decomposed into

with

where $\theta_{i}$ represents the angle of the Gaussian function’s major axis, and $\sigma^{\text{maj}}_{i}$ and $\sigma^{\text{min}}_{i}$ represent its extent in major and minor axis, respectively.

The final loss function for simultaneously regressing $N$ heatmaps is defined as

Same as in Eq. (3), the first term in Eq. (7) minimizes the differences between the predicted heatmaps $\hat{h}_{i}$ and the target heatmaps $h_{i}$ for all landmarks $L_{i}$, while we now not only learn the isotropic $\sigma_{i}$ but all covariance parameters of the Gaussian functions ($\theta_{i}$, $\sigma^{\text{maj}}_{i}$, and $\sigma^{\text{min}}_{i}$). Same as before, the second term of Eq. (3) is used to prefer smaller Gaussian extents by penalizing $\sigma^{\text{maj}}_{i}\sigma^{\text{min}}_{i}$ with factor $\alpha$. Note that while the mean parameters of the Gaussian functions ($\mathbf{x}_{i}$) are set to the groundtruth landmark locations and differ for each training image, the covariance parameters of the target heatmaps $h_{i}$ are learned for the whole dataset.

Analyzing Eq. (7) in more detail, not only the predicted heatmaps $\hat{h}_{i}$ aim to be close to the target heatmaps $h_{i}$ but also the target heatmaps $h_{i}$ aim to be close to the predicted heatmaps $\hat{h}_{i}$. For each landmark $L_{i}$, $\hat{h}_{i}$ and $h_{i}$ receive feedback from each other during training (see Fig. 1). While the network parameters $\mathbf{w}$ are updated to better model the shape of the target heatmap $\hat{h}_{i}$, at the same time, the covariance $\mathbf{\Sigma}_{i}$ parameters of each target heatmap $h_{i}$ ($\theta_{i}$, $\sigma^{\text{maj}}_{i}$, and $\sigma^{\text{min}}_{i}$) are updated to better model the shape of the predicted heatmap $\hat{h}_{i}$. Note that, while $\mathbf{\Sigma}_{i}$ is learned from a single annotation per image, the CNN sees one annotation for each image in the training set. This allows implicitly modeling the ambiguities of landmark annotations for the whole training dataset representing a dataset-based annotation uncertainty.

As the output of heatmap regression networks are $N$ heatmaps that correspond to the $N$ landmarks $L_{i}$ rather than their coordinates, each landmark’s coordinate $\hat{\mathbf{x}}_{i}$ needs to be obtained from the predicted heatmaps $\hat{h}_{i}$. Previous work often uses the coordinate of the maximum response or the center of mass of the predicted heatmaps as the landmark’s coordinate. However, this represents only the most probable landmark location, but not the whole distribution of possible landmark locations.

As the target heatmaps are modeled as Gaussian functions, the predicted heatmaps are also expected to be Gaussian functions, see Eq. (4). Hence, the distribution of possible landmark locations is captured by fitting Gaussian functions to the predicted heatmaps $\hat{h}$. We use a robust least squares curve fitting method (Branch et al., 1999) to fit the Gaussian function (see Eq. (4)) to the predicted heatmaps $\hat{h}_{i}$ initialized to the respective maximum and obtain the fitted Gaussian parameters $\hat{\mathbf{x}}_{i}$ and $\hat{\mathbf{\Sigma}}_{i}$ with $i=\{1,...,N\}$. As the fitted Gaussian function models the distribution of possible locations, the parameters represent a directional sample-based annotation uncertainty of the landmark location during inference.

The proposed method is evaluated on two publicly available datasets for anatomical landmark localization. The first dataset consists of 895 radiographs of left hands of the Digital Hand Atlas Database System (Gertych et al., 2007; Zhang et al., 2009) with 37 annotated landmarks per image. Each image has been annotated by one of three experts. The second dataset consists of 400 lateral cephalograms with 19 annotated landmarks per image. This dataset has been used for the ISBI 2015 Cephalometric X-ray Image Analysis Challenge (Wang et al., 2016), in which each image has been annotated by both a senior and a junior radiologist.

For evaluating the landmark localization performance, the point-to-point error is used for the landmarks $i$ in image $j$, which is defined as the Euclidean distance between the target coordinate $\mathbf{x}^{(j)}_{i}$ and the predicted coordinate $\hat{\mathbf{x}}^{(j)}_{i}$, i.e., $\text{PE}_{i}^{(j)}=\lVert\mathbf{x}^{(j)}_{i}-\hat{\mathbf{x}}^{(j)}_{i}\rVert_{2}$. This allows calculation of sets of point-to-point errors, e.g., for all $M$ images of the dataset over all $N$ landmarks, over which the mean and standard deviation can be calculated. Additionally, we report the success detection rate, which is defined as the percentage of predicted landmarks that are within a certain point-to-point error radius $r$ for all images, i.e., $\text{SDR}_{r}=\frac{100}{MN}\big{|}\{\,(j,i)\;|\;\text{PE}_{i}^{(j)}\leq r\}\big{|}$.

While in most experiments only a single groundtruth coordinate per image and landmark is available, the groundtruth location in our inter-observer experiment is not only encoded as a single coordinate, but as a distribution. Such groundtruth distributions are compared to the distributions predicted by our method by analyzing the parameters of the Gaussian functions. As a measure of size of the distributions, we use the area encoded by the product of the extent of the major axis and minor axis, i.e., $\sigma^{\text{maj}}\cdot\sigma^{\text{min}}$. As a measure of anisotropy, we use their ratio, i.e., $\sigma^{\text{maj}}:\sigma^{\text{min}}$. Additionally, to compare the predicted and groundtruth orientation for landmarks with a strong anisotropic distribution, the orientation angle $\theta$ with respect to the x-axis is reported.

In Tables 1 and 2, landmark prediction accuracy measured in point error (PE) and robustness measured in the success detection rate per radius ($\text{SDR}_{r}$) of our method is compared to the state-of-the-art on both the hand and cephalogram dataset, respectively. The reported values for the state-of-the-art methods are taken from the respective publications (Ibragimov et al., 2014; Lindner et al., 2015; Urschler et al., 2018; Payer et al., 2019; Zhong et al., 2019; Chen et al., 2019).

On a three-fold cross-validation (CV) of the dataset of left hand radiographs, our proposed method that uses an anisotropic Gaussian function to model both the target heatmap and the distribution of the heatmap prediction ($\mathbf{\Sigma}$-target,$\hat{\mathbf{\Sigma}}$-fit) outperforms the previously published methods (Lindner et al., 2015; Urschler et al., 2018; Payer et al., 2019). Also, our method is compared to the winning methods of the ISBI Cephalometric X-ray Image Analysis Challenges in 2014 (Ibragimov et al., 2014) and 2015 (Lindner et al., 2015) as well as recent state-of-the-art methods on this dataset (Zhong et al., 2019; Chen et al., 2019). On both test sets Test1 and Test2, it can be observed that the methods using CNNs achieve similar results, while our method is the most accurate with the smallest PE. In terms of SDR, our method performs in-line with other methods, while Chen et al. (2019) perform slightly better in Test1 and Zhong et al. (2019) in Test2, respectively. Nevertheless, the results on Test1 and Test2 should be treated with caution because of the systematic shift of some landmarks (e.g. landmark 8 and 16) as compared to the annotations of the training data. The systematic shift has already been reported in  (Lindner et al., 2016; Zhong et al., 2019) as well as in our MICCAI-UNSURE 2020 workshop paper (Payer et al., 2020). To mitigate the systematic shift, we follow previous work (Lindner et al., 2016; Zhong et al., 2019) and perform a four-fold cross-validation experiment trained and evaluated on the junior annotations only (CV jun.). Here, similar to the hand radiographs, it can be observed that our method that uses an anisotropic Gaussian function to model both the target heatmap and the distribution of the heatmap prediction ($\mathbf{\Sigma}$-target,$\hat{\mathbf{\Sigma}}$-fit) leads to the best results, outperforming all other methods.

In this work, we proposed to use parameters of the target Gaussian heatmaps to model a dataset-based annotation uncertainty (see Sec. 3.1) and parameters of a fitted Gaussian function to the predicted heatmaps to model a sample-based annotation uncertainty (see Sec. 3.2). In this section, it is shown that both sets of Gaussian parameters represent valid uncertainty measures by comparing the sizes of the Gaussians with the corresponding localization error for each landmark individually (see Fig. 3). For both the hand and cephalogram dataset, the size of the learned target Gaussians as well as the size of the fitted Gaussians are calculated as the product of the variances ($\sigma^{\text{maj}}\cdot\sigma^{\text{min}}$) averaged over the cross-validations or over all test images, respectively. It can be observed that both the target and fitted heatmap parameters correlate with the PE.

Our proposed method is also able to model directional uncertainties, which is visualized in Fig. 4 and 5 for selected landmarks. On the left side of the figures, the target and predicted heatmaps are superimposed with ellipses representing the target and fitted Gaussian parameters, respectively. Both target and predicted heatmaps are smaller for landmarks on distinct corners as compared to landmarks on edges. For landmarks on edges it can be observed that the heatmaps are oriented in direction of the underlying edges. Thus, for these landmarks the main sources for ambiguity are along the edges. In Fig. 4 and 5 (right) it is shown that both the learned target heatmap parameters $\mathbf{\Sigma}_{i}$ and the fitted heatmap parameters $\hat{\mathbf{\Sigma}}_{i}$ correspond in orientation and size to the expected prediction-groundtruth offsets $\mathbf{x}_{i}-\hat{\mathbf{x}}_{i}$. For each landmark individually, the parameters of both target and predicted heatmaps are represented as ellipses centered at the origin, and plotted together for all images. To visualize the prediction-groundtruth offsets, individual predictions for all images are shown as points relative to their respective groundtruth positions, centered at the origin. From Fig. 4 and 5 it can be observed that the orientation and the size of both predicted and target heatmap parameters correspond to the prediction-groundtruth offsets. Overall, these results confirm that the target and the fitted parameters of the Gaussian functions are valid uncertainty measures.

To reliably assess inter-observer variability, neither the hand dataset nor the cephalogram dataset can be used, as they were annotated by only one or two annotators per image, respectively. Therefore, we randomly chose 100 images of the cephalogram dataset and employed nine annotators experienced in medical image analysis to label five characteristic landmarks on each image. Combined with the junior and senior radiologist annotations already existing in the cephalogram dataset, this leads to a total of 11 annotations per image and landmark.

We trained a model of our proposed method ($\mathbf{\Sigma}$-target,$\hat{\mathbf{\Sigma}}$-fit) on the mean coordinate according to the 11 annotations per image. The annotation and predicted distributions of the selected points are compared in Fig. 6. Similarly to how the prediction-groundtruth offsets are presented in Fig. 4, the annotation offsets for each annotator and image are visualized as points relative to the respective mean annotation of all 11 annotators, centered at the origin (Fig. 6, first column). To model the annotation distribution, a Gaussian distribution is fitted to the 11 annotations for each landmark and image individually. All annotation distributions as well as their mean are visualized as ellipses (Fig. 6, second column). The parameters of both predicted heatmaps (Fig. 6, third column) and target heatmaps (Fig. 6, last column) are represented as ellipses centered at the origin, plotted for all images. By visually comparing the size and the orientation of both predicted and target heatmap parameters with the annotation distributions, a correspondence can clearly be seen.

When comparing the annotation offsets for different landmarks in Fig. 6, it can be observed that landmark 4 on the tip of the incisor has the best agreement among the annotators, due to its well defined location (see Fig. 2). Although landmark 5 also lies on an anatomically well defined location, due to low contrast in this region, its annotation variation is higher. Landmarks 2 and 3 were the most problematic ones to annotate. Landmark 2 is challenging, because its position has to be estimated by the annotators in between the left and right mandible, when they are not aligned and therefore, are not seen as a single edge in the lateral cephalograms. Landmark 3 lies on an edge of a bone within the skull that is challenging to recognize and becomes ambiguous if the head is tilted. For landmark 1 that lies on the tip of the chin, it can be seen that the annotators correctly identified the edge of the chin, however, identifying the exact position along the edge of the chin was challenging. Thus, the annotation distribution for this landmark has a clear dominant direction.

In all the previous works based on heatmap regression, the target heatmaps are modeled as isotropic functions. However, from the annotation distributions in Fig. 6, second column, clear anisotropic orientations of the inter-observer variabilities of the landmarks located at smooth edges, i.e., landmark 1 and 2 are observable. Thus, to estimate the inter-observer variability, in this work, we used anisotropic Gaussian functions for target heatmaps, as well as to model the predicted distribution. This enables our method to successfully model the inter-observer variability as can be seen when comparing both the predicted Gaussian functions in the third and the learned target Gaussian functions in the last column to the annotation distributions in the second column of Fig. 6.

To quantitatively validate this correspondence between the parameters of the Gaussian functions of our proposed method and the annotation distribution, we present the ratio and product of the extent of the major axis $\sigma^{\text{maj}}$ and minor axis $\sigma^{\text{min}}$ as well as the orientation angle $\theta$ of the Gaussian functions in Table 3. The product of the extents of the major axis $\sigma^{\text{maj}}$ and minor axis $\sigma^{\text{min}}$ shows that the parameters of both target and average fitted Gaussians have the tendencies of underestimating the size of the annotation uncertainty for landmarks with a larger inter-observer variability (landmark 2, 3, 5). Nevertheless, for the landmarks with lower inter-observer variability (landmark 1, 4), the size of the Gaussian functions is in range with the size of the annotation uncertainty. From the ratio of the extent of the major axis $\sigma^{\text{maj}}$ and minor axis $\sigma^{\text{min}}$, it can be observed that our method is not only able to model isotropic annotation distributions (landmark 3, 4, 5) but also anisotropic ones (landmark 1, 2). It is interesting to see that both ratios of the target and the average fitted Gaussian parameters are similar to the ratio of the annotation distributions. Furthermore, the rotation angles $\theta$ of the Gaussian functions also approximate the rotation angle of the annotation distribution for the anisotropic landmarks well.

Qualitative results of the landmark predictions as well as the predicted distributions for each landmark are presented in Fig. 7 for selected images. The visualized points are the 11 inter-observer annotations (dark red), their mean position (red) and the landmark predictions of our model (blue). The red ellipses represent the distributions fitted to the 11 annotations. The sample-based uncertainties and the dataset-based uncertainties are represented as blue and green ellipses, respectively. For the depicted images, it can be seen that most landmark predictions closely match to the mean of the 11 annotations, while the predicted uncertainties are similar to the annotation distributions. Discrepancy between the inter-observer annotations and uncertainties predicted by our method were primarily caused by landmarks that where difficult to recognize in a given image (landmark 2, bottom left and landmark 4, bottom middle) or by landmarks with high ambiguity (e.g., landmark 3, bottom middle and landmark 2 bottom right). Nevertheless, the predicted sample-based uncertainties are similar to the annotation distributions even when the dataset-based uncertainty of a landmark is too restrictive as, e.g., for landmark 2 (top middle), where both mandibles are visible due to the head being tilted.

Different strategies of using the Gaussian function to model both the target heatmap and the distribution of the heatmap prediction are investigated. We distinguish methods based on whether the target heatmap used during training is defined as a Gaussian function with a fixed isotropic variance ($\sigma$=3), learned isotropic variance ($\sigma$-target) or learned anisotropic covariance ($\mathbf{\Sigma}$-target). For the fixed isotropic variance, the value $\sigma=3.0$ was chosen, as it showed to have a good trade-off between PE and $\text{SDR}_{r}$. As already defined, our main method uses the Gaussian function with anisotropic covariance during both training and inference, which we named $\mathbf{\Sigma}$-target,$\hat{\mathbf{\Sigma}}$-fit. For comparison with state-of-the-art methods estimating uncertainty, we implemented two approaches using MCD as described in Sec. 4. The cephalogram dataset with our additional annotations was used for method comparison as well as the same visualizations and measures as in the previous Sec. 5.3. Results of these comparisons for two selected landmarks with the highest anisotropy are shown in Fig. 8, while the quantitative results for all five landmarks are shown in Table 4.

When comparing the localization accuracy in terms of PE, all different heatmap target and fitting strategies perform similarly well as can be seen in Table 4. Similarly as in Sec. 5.3, it can be seen that both target and predicted heatmaps have the tendency of underestimating the annotation distribution size $\sigma^{\text{maj}}\cdot\sigma^{\text{min}}$ for landmarks having a large inter-observer variability. When looking at the distributions modeled by the target heatmaps in Fig. 8, it can be observed that strategies using either fixed ($\sigma=3$) or learned isotropic $\sigma$ values ($\sigma$-target) have no capacity of modeling the annotation distribution. Interestingly, when fitting anisotropic Gaussian functions to the heatmaps predicted by both of these strategies, the orientation angles $\theta$ of the inter-observer variability are successfully estimated. However, the ratio of $\sigma^{\text{maj}}:\sigma^{\text{min}}$ in Table 4 shows that our final model, which uses anisotropic Gaussian functions for both target heatmaps and prediction fitting ($\mathbf{\Sigma}$-target,$\hat{\mathbf{\Sigma}}$-fit), better captures the anisotropy of the annotator distributions.

Our method is also compared with the state-of-the-art method MCD for estimating uncertainties adapted to the landmark localization task, where the uncertainty is either computed as the variance of a set of stochastic forward predictions (MCD-$\hat{\mathbf{x}}$,max) or estimated by fitting a Gaussian to an accumulated set of stochastic heatmaps (MCD-$\hat{h}$,$\hat{\mathbf{\Sigma}}$-fit). Since fixed target heatmaps were used to be consistent with related work, only the variances respectively the parameters of the fitted Gaussians are compared with the annotation distributions. As can be seen in Table 4, there is no significant difference when comparing MCD to any of our strategies in terms of the PE. However, MCD-$\hat{\mathbf{x}}$,max is not capable of estimating the size of the annotation distributions $\sigma^{\text{maj}}\cdot\sigma^{\text{min}}$, due to the predicted coordinates being very similar for all passes which leads to an underestimation. Moreover, the ratio of $\sigma^{\text{maj}}:\sigma^{\text{min}}$ is also unreliable for this method. The method MCD-$\hat{h}$,$\hat{\mathbf{\Sigma}}$-fit has a similar performance as our strategy $\sigma$-target,$\hat{\mathbf{\Sigma}}$-fit, which was to be expected as the uncertainty is also estimated by fitting a Gaussian.

In clinical practice utilizing computer aided diagnosis, often only the most probable landmark prediction is used, while the uncertainty distribution of the prediction is ignored leading to overconfidence in calculated clinical measurements. To show the importance of considering the annotation uncertainty of the dataset used for training the machine learning model, we evaluate the influence of the predicted landmark localization uncertainty to the classification accuracy of clinical measurements for orthodontic analysis and treatment planning. For this experiment, the same eight measurements that classify anatomical abnormalities as described in Wang et al. (2016) and listed in Fig. 9 are used.

To show the correlation of the predicted landmark localization uncertainty with the classification accuracy, we first calculated how the classification prediction is affected by the predicted landmark distributions. Thus, for each patient and measurement, 10,000 points are sampled from each landmark distribution predicted by our $\mathbf{\Sigma}$-target,$\hat{\mathbf{\Sigma}}$-fit method trained using the junior annotations to calculate 10,000 classification outcomes in the four-fold cross-validation setup (CV jun.) that we used before in Sec. 5.1. For each patient and measurement, the entropy of the class probabilities is used to obtain the classification uncertainty, while the classification accuracy is calculated as the agreement of the most probable class prediction with the junior annotation.

To verify that such a calculated entropy is a valid measure of uncertainty, we compare the classification uncertainty with the classification accuracy, for each clinical measurement independently. Fig. 9 visualizes the average accuracy that is calculated by increasing the number of considered patients, ordered from lower to higher uncertainty. For all clinical measurements, the blue line in Fig. 9 shows that the average classification accuracy is decreasing when the number of images with higher classification uncertainty is increasing. Thus, our experiment shows that the predicted landmark localization uncertainty can also be used to model classification uncertainty influencing the classification accuracy.

Additionally, the dashed line in Fig. 9 shows the classification accuracy of the senior as compared to the junior annotator, for the same order of images as used to generate the blue line. When comparing the dashed with the blue line, it can be seen that the classification agreement of the senior to the junior is in general lower as compared to our method, which can be explained by our model being trained on the junior annotation. However, it is more interesting to see that by ordering the images according to the uncertainty of our model, the accuracy of the senior’s classification compared to the junior’s is also decreasing. An interesting extension would be to estimate the optimal heatmap size directly by maximizing the classification accuracy, however, this would go beyond the scope of this work focusing on landmark localization. Nevertheless, the experiment in this subsection further confirms that our method was able to identify samples with higher ambiguities in annotation.