Relational reasoning network for anatomical landmarking

Abstract. Purpose We perform anatomical landmarking for craniomaxillofacial (CMF) bones without explicitly segmenting them. Toward this, we propose a simple, yet efficient, deep network architecture, called relational reasoning network (RRN), to accurately learn the local and the global relations among the landmarks in CMF bones; specifically, mandible, maxilla, and nasal bones. Approach The proposed RRN works in an end-to-end manner, utilizing learned relations of the landmarks based on dense-block units. For a given few landmarks as input, RRN treats the landmarking process similar to a data imputation problem where predicted landmarks are considered missing. Results We applied RRN to cone-beam computed tomography scans obtained from 250 patients. With a fourfold cross-validation technique, we obtained an average root mean squared error of <2  mm per landmark. Our proposed RRN has revealed unique relationships among the landmarks that help us in inferring informativeness of the landmark points. The proposed system identifies the missing landmark locations accurately even when severe pathology or deformations are present in the bones. Conclusions Accurately identifying anatomical landmarks is a crucial step in deformation analysis and surgical planning for CMF surgeries. Achieving this goal without the need for explicit bone segmentation addresses a major limitation of segmentation-based approaches, where segmentation failure (as often is the case in bones with severe pathology or deformation) could easily lead to incorrect landmarking. To the best of our knowledge, this is the first-of-its-kind algorithm finding anatomical relations of the objects using deep learning.


Introduction
In the United States alone, more than 17 million patients suffer from developmental deformities of the jaw, face, and skull region due to trauma, deformities from tumor ablation, or congenital birth defects. 1 The number of patients who require orthodontic treatment is far beyond this number. An accurate anatomical landmarking on medical scans (mostly it is volumetric com- puted tomography-CT-scans) is a crucial step in the deformation analysis and surgical planning of the craniomaxillofacial (CMF) bones. This, if done correctly and efficiently, would afford precise image-based surgical planning for patients. This is even more significant since such deformities are known to vary from patient to patient and hence need careful delineation. However, manual landmarking in volumetric CT scans is a tedious process and prone to inter-operator variability. There are considerable efforts towards developing a fully-automated and accurate software for anatomical landmarking based on bone segmentation from CT scans. [2][3][4] Despite this clinical need, very little progress has been made especially for bones with a high level of congenital and developmental deformations (approximately 5% of the CMF deformities). Deep learning based approaches have become the standard choice for pixel-wise medical-image segmentation applications due to their high efficacy. 2,5,6 However, it is difficult to generalize segmentation especially when there is a high degree of deformation or pathology, 7 which is the case for treating CMF conditions. intervention (left) and high variability in the bone (right), and causing segmentation algorithms to fail (leakage or under-segmentation). Current state-of-the-art landmarking algorithms are mostly dependent on bone segmentation results, since locating landmarks could become easier once their parent anatomy (the bones they belong to) is precisely known. 7 Figure 2 demonstrates mandible and maxilla/nasal bone anatomies along with the landmarks associated with those bones. If the underlying segmentation is poor, it is highly likely to have a large landmark localization error, directly affecting the quantification process (which could include severity measurement, surgical modeling, and treatment planing).
We hypothesize that if an explicit segmentation can be avoided for extremely challenging cases, landmark localization errors can be minimized. This will also lead to a widespread use of landmarking procedure in surgical planning and precision medicine. Since, CMF bones are found in the same anatomical space even when there is deformity or pathology. Hence, the overall global relationship of the anatomical landmarks should still be preserved despite severe localized changes.
Based on this rationale, we claim that utilizing local and global relations of the landmarks can help automatic landmarking without the extreme need for segmentation.

Landmarking
Anatomical landmark localization approaches can broadly be categorized into three main groups: registration-based (atlas-based), 8 knowledge-based, 9,10 and learning-based. 7,11 Integration of shape and appearance increases the accuracy of the registration-based approaches. However, image registration is still an ill-posed problem, and when there are variations such as age (pediatrics vs. adults), missing teeth (very common in certain age groups), missing bone or bone parts, severe pathology (congenital or trauma), and imaging artifacts, the performance can be quite poor. 3,12,13 The same concerns apply to segmentation based approaches too.
Gupta et al. 10 developed a knowledge-based algorithm to identify 20 anatomical landmarks on cone-beam CT (CBCT) scans. Despite their promising results, a seed must be selected by using 3D template registration on the inferior-anterior region where fractures are most commonly found. An error in the seed localization may easily lead to a sub-optimal outcome in such approaches. Zhang et al. 14 developed a regression forest-based landmark detector to localize CMF landmarks on the CBCT scans. To address the spatial coherence of landmarks, image segmentation was used as a helper. The authors obtained a mean digitization error less than 2 mm for 15 CMF landmarks. The following year, to reduce the mean digitization error further, Zhang et al. 2 proposed a deep learning based joint CMF bone segmentation and landmarking strategy. A context guided multi-task fully convolutional neural (FCN) network was employed along with 3D displacement maps to perceive the spatial locations of the landmarks. A segmentation accuracy of 93.27 ± 0.97% and a mean digitization error of less than 1.5 mm for identifying 15 CMF landmarks was achieved. Further, a joint segmentation and landmark digitization framework was proposed, where two stages of FCN were cascaded to perform bone segmentation and landmark localization. 7 The major disadvantage of this (one of the state-of-the-arts) method was the memory constraint introduced by the redundant information in the 3D displacement maps such that only a limited number of the landmarks can be learned using this approach. Since the proposed strategy is based on joint segmentation and landmarking, it naturally shares other disadvantages of the segmentation based methods: frequent failures for very challenging cases. The landmark localization problem was solved using an object detection method, where region proposals were used to identify landmark locations and a coarse to fine method was used to achieve landmark localization. 15 It must be noted, that the method does not use the relationships between the anatomical landmarks in the CMF bones.
Recently, we integrated the manifold information (geodesic) in a deep learning architecture to improve robustness of the segmentation based strategies for landmarking, 5 and obtained promising results, significantly better than the state-of-the-art methods. We also noticed that there is still a room to improve landmarking process especially when pathology or bone deformation is severe.
To fill this research gap, in this study, we take a radically different approach by learning landmark relationships without segmenting bones. We hypothesize that the inherent relation of the landmarks in the CMF region can be learned by a relational reasoning algorithm based on deep learning.
Although our proposed algorithm stems from this unique need of anatomical landmarking, the core idea of this work is inspired from the recent studies in artificial intelligence (AI), particularly in robotics and physical interactions of human/robots with their environments, as described in the following with further details.

Relational Reasoning
The ability to learn relationship and infer reasons between entities and their properties is a central component of the AI field, however it has been proven to be very difficult to learn through neural networks until recently. 16,17 In 2009, Scarselli et al. 18 introduced a graph neural network (GNN) by extending the neural network models to process graph data which encoded relationship information of the objects under investigation. Li et al. 19 proposed a machine learning model based on gated recurrent units (GRUs) to learn the distributed vector representations from heap graphs. Despite this increase in use and promising nature of the GNN architectures, 20 there is a limited understanding for their representational properties, which is often a necessity in medical AI applications for their adoption in clinics.
Recently, DeepMind team(s) published four important studies on relational reasoning and explored how objects in complex systems can interact with each other. 16, 21-23 Battaglia et al. 21 introduced interaction networks to reason about the objects and the relations in the complex environments. The authors proposed a simple, yet accurate system to reason about n-body problems, rigid-body collision, and non-rigid dynamics. The proposed system can predict the dynamics in the next step with an order of magnitude lower error and higher accuracy. Raposa and Santoro et al. 16 introduced a Relational Network (RN) to learn the object relations from a scene descrip-tion, hypothesising that a typical scene contains salient objects which are typically related to each other by their underlying causes and semantics. Following this study, Santoro and Raposa et al. 22 presented another relational reasoning architecture for tasks such as visual question-answering, text-based question-answering, and dynamic physical systems where the proposed model obtained most answers correctly. Lastly, Battaglia et al. 23 studied the relational inductive biases to learn the relations of the entities and presented the graph networks. These four studies show promising approaches to understanding the challenge of relational reasoning. To the best of our knowledge, such advanced reasoning algorithms have neither been developed for nor applied to the medical imaging applications yet. It must be noted that medical AI applications require fundamentally different reasoning paradigms than conventional computer vision and robotics fields have 24 (e.g., salient objects definitions). To address this gap, in this study we focus on the anatomy-anatomy and anatomy-pathology relationships in an implicit manner.

Summary of our contributions
• The proposed method is the first in the literature to successfully apply spatial reasoning of the anatomical landmarks for accurate and robust landmarking using deep learning.
• Many anatomical landmarking methods, including our previous works, 5, 10, 25 use bone segmentation as a guidance for finding the location of the landmarks on the surface of a bone.
The major limitation imposed by such an approach stems from the fact that it is not always possible to have an accurate segmentation. Our proposed RRN system addresses this problem by enabling accurate prediction of anatomical landmarks without employing explicit object segmentation.
• Since efficiency is a significant barrier for many medical AI applications, we explore new deep learning architecture designs for a better efficacy in the system performance. For this purpose, we utilize variational dropout 26 and targeted dropout 27 in our implementation for faster and more robust convergence of the landmarking procedure (∼ 5 times faster than baselines).
• Our data set includes highly variable bone deformities along with other challenges of the CBCT scans with a larger number of scans (as compared to baselines). Hence, the proposed algorithm is considered robust and identifies anatomical landmarks accurately under varying conditions (Table 3). In our experiments, we find landmarks pertaining to mandible, maxilla, and nasal bones ( Figure 2). The rest of this paper is organized as follows: we introduce our novel methodology and its details in Section 2. In Section 3, we present experiments and results and then we conclude the paper with discussing strengths and limitations of our study in Section 4.

Overview and preliminaries
The most frequently deformed or injured CMF bone is the lower jaw bone, mandible, which is the only mobile CMF bone. 28 In our previous study, 5 we developed a framework to segment mandible from CBCT scans and identify the mandibular landmarks in a fully-automated way. Herein, we focus on anatomical landmarking without the need for explicit segmentation, and extend the learned landmarks into other bones (maxilla and nasal). Overall, we seek answers for the following important questions: • Q1: Can we automatically identify all anatomical landmarks of a bone if some of the landmarks are missing? If so, what is the least effort for performing this procedure? How many landmarks are necessary and which landmarks are more informative to perform this whole procedure?
• Q2: Can we identify anatomical landmarks of nasal and maxilla bones if we only know locations of a few landmarks in the mandible and the rest is missing? Do relations of landmarks hold true even when they belong to different anatomical structures (manifold)?
In this study, we explore inherent relations among anatomical landmarks at the local and global levels in order to explore availability of structured data samples helping anatomical landmark localization. Inferred from the morphological integration of the CMF bones, we claim that landmarks of the same bone should carry common properties of the bone so that one landmark should give clues about the positions of the other landmarks with respect to a common reference. This reference is often chosen as segmentation of the bone to enhance information flow, but in our study we leverage this reference point from the whole segmented bone into a reference landmark point. Throughout the text, we use the following definitions: Definition 1: A landmark is an anatomically distinct point, helping clinicians to make reliable measurements related to a condition, diagnosis, modeling a surgical procedure, or creating a treatment plan.

Definition 2:
A relation is defined as a geometric property between landmarks. Relations might include the following geometric features: size, distance, shape, and other implicit structural information. In this study, we focus on pairwise relations between landmarks as a starting point. Once relationship among landmarks is learned effectively, we can use this relationship to identify the missing landmarks on the same or different CMF bones without the need for a precise segmentation. Towards this goal, we propose to learn a relationship of the anatomical landmarks in two stages (illustrated in Figure 3) based on Relational Units (RUs). The first stage which is shown as the function g, which learns the pairwise (local) relations. The second stage which is shown as a function f , which combines pairwise relations (g) of the first stage into a global relation based on RUs.

Relational Reasoning Architecture
Anatomical landmarking has been an active research topic for several years in the medical imaging field. However, how to build a reliable/universal relationship between landmarks for a given clini-  noted that this combination is employed through a joint loss function and an RU to infer an average relation information. In other words, for each individual landmark, the combined relationship vector is assigned a secondary learning function through a single RU.
The RU is the core component of the RRN architecture. Each RU is designed in an end-to-end fashion; hence, they are differentiable. For n landmarks in the input domain, the proposed RRN architecture learns n × (n − 1) pairwise and n combined relations (global) with a total of n 2 RUs.
Therefore, depending on the number of input domain landmarks, RRN can be either shallow or dense. Let L input andL indicate vectors of input and output anatomical landmarks, respectively.
Then, two stages of the RRN of the input domain landmarks L input can be defined as:  Table 2). Since all pairwise relations are leveraged under G θ i and f φ with averaging operation, we can conclude that RRN is invariant to the order of input landmarks (i.e., permutation-invariant).

Loss Function
The natural choice for the loss function is the mean squared error (MSE) because it is a differentiable distance metric measuring how well landmarks are localized/matched, and it allows output of the proposed network to be real-valued functions of the input landmarks. For n input landmarks and m target landmarks, MSE simply penalizes large distances between the landmarks as follows: where o k are target domain landmarks (o k ∈L).

Variational Dropout
Dropout is an important regularizer employed to prevent overfitting at a cost of 2 − 3 times (on average) increase in training time. 32 For efficiency reasons, speeding up dropout is critical and it can be achieved by a variational Bayesian inference on the model parameters. 26 Given a training input dataset X = {x 1 , x 2 , .., x N } and the corresponding output dataset Y = {y 1 , y 2 , .., y N }, the goal in RRN is to learn the parameters ω such that y = F ω (x). In the Bayesian approach, given the input and output datasets X, Y , we seek for the posterior distribution p(ω|X, Y ), by which we can predict output y * for a new input point x * by solving the integral: 33 In practice, this computation involves intractable integrals. 26 To obtain the posterior distributions, a Gaussian prior distribution N (0, I) is placed over the network weights 33 which leads to a much faster convergence. 26

Targeted Dropout
Alternatively, we also propose to use targeted dropout for better convergence. 27 Given a neural network parameterized by Θ, the goal is to find the optimal parameters W Θ (.) such that the loss Loss(W Θ ) is minimized. For efficiency and generalization reasons, |W Θ | ≤ k, only k weights of highest magnitude in the network are employed. In this regard, deterministic approach is to drop the lowest |W Θ | − k weights. In targeted dropout, using a target rate γ and a drop out rate α, first a target set is generated with the lowest weights with the target rate γ. Next, weights are dropped out in an stochastic manner from the target set at a certain dropout rate α.    (Table 1). Finally, a 19-dimensional feature vector is considered to be an input to local relationship function g. For a well-defined reference landmark, we used Menton (Me) as the origin of the Mandible region.

Data Description
Anonymized CBCT scans of 250 patients (142 female and 108 male, mean age = 23.6 years, standard deviation = 9.34 years) were included in our analysis through an IRB-approved protocol.
The data set includes both pediatric and adult patients with craniofacial congenital birth defects, developmental growth anomalies, trauma to the CMF, and surgical interventions. CB MercuRay In addition, following image-based variations exist in the data set: aliasing artifacts due to braces, metal alloy surgical implants (screws and plates), dental fillings, and missing bones or teeth. 5 The data was annotated independently by three expert interpreters, one from the NIH team, and two from UCF team. Inter-observer agreement values were computed as approximately 3 pixels.
Experts used freely available 3D Slicer software for the annotations. 5

Data Augmentation
Our data set includes fully-annotated mandibular, maxillary, and nasal bones' landmarks. Due to insufficiency of 250 samples for a deep-learning algorithm to run, we applied data-augmentation approach. In our study, the common usage of random scaling or rotations for data-augmentation were not found to be useful for new landmark data generation because such transformations would not generate new relations different from the original ones. Instead, we used random interpolation similar to active shape model's landmarks. 34 Briefly, we interpolated 2 (or 3) randomly selected scans with randomly computed weight. We merged the relation information at different scans to a new relation. We also added random noise to each landmark with a maximum in the range of ±5 pixels, defined empirically based on the resolution of the images as well as the observed highdeformity of the bones. We generated a dataset with approximately 100K landmarks, which is empirically evaluated as a sufficiently large dataset.

Evaluation Methods
We used root-mean squared error (RMSE) in the anatomical space (in mm) to evaluate the goodness of the landmarking. Lower RMSE indicates successful landmarking process. For statistical significance comparisons of different methods and their variants, we used a P-value of 0.05 as a cut-off threshold to define significance and applied t-tests where applicable.

Input Landmark Configurations
In our experiments, there were three groups of landmarks (See Figure 2)  • L 2 = {Ans, A, P r, P ns, }, In each experiment, as detailed in Table 2, we designed a specific input set L input where in Section 2: • What configuration of the input landmarks can capture the manifold of bones so that other landmarks can be localized successfully?
• What is the minimum number and configuration of the input landmarks for successful identification of other landmarks?

Training
The MLP RU was composed of 3 fully-connected layers, 2 batch normalizations and 2 ReLUs

Experiments and Results
We ran a set of experiments to evaluate the performance of the proposed system using a 4-fold cross-validation. We summarized the experimental configurations in Table 2, error rates in Table 3, and corresponding renderings in Figure 6. The method achieving the minimum error for a corresponding landmark is colored the same as the corresponding landmark at Table 3. As shown by results, the minimum number of the input landmarks required for successful identification of other landmarks is determined as 3.
Among two different RU architectures, DB architecture was evaluated to be more robust and fast to converge as compared to the MLP architecture. To be self-complete, we provided the MLP experimental configuration performances only for the 5-landmark experiment (See Table 3).
In the first experiment (Table 2-Experiment 1), to have an understanding of the performance of the RRN, we used the landmark grouping sparsely-spaced and closely-spaced as proposed in Torosdagli et al. 5 We named our first configuration as "5-landmarks" where closely-spaced, max-  Table 2). In the 5-landmarks RRN architecture, there were totally 25 RUs. In the second experiment (Table 2-Experiment 2), we explored the impact of a configuration with less number of input mandibular landmarks on the learning performance. Compared to the 5 sparsely-spaced input landmarks, we learned the relation of 3 landmarks, M e, Cd L and Cd R , and predicted the closelyspaced landmark locations (as in the 5-landmarks experiment) plus superior-anterior landmarks Cor L and Cor R and maxillary and nasal bones' landmark locations. The network was composed of 9 RUs. The training was relatively fast compared to the 5-landmarks configuration due to small number of RUs. We named this method as "3-Landmarks Regular".
After observing statistically similar accuracy compared to the 5-landmarks method for the closely-spaced landmarks (P > 0.005), and high error rates at the superior-anterior landmarks Cor L and Cor R , we setup a new experiment which we named "3-Landmarks Cross" (Table 2-Experiment 3). In this configuration, the third experiment, we used 1 superior-posterior and 1 superior-anterior landmarks on the right and left sides, respectively. This network was similar to 3-landmarks regular one in terms of number of RUs used.
In the fourth experiment (Table 2-Experiment 4), we evaluated the performance of the system in learning the closely-spaced mandibular landmarks (Gn, P g, B, Id) and the maxillary landmarks (AN S, A, P r, P N S) using the relation information of the sparsely-spaced and the nasal-bones landmarks which is named as "6-landmarks". There are a total of 36 RUs in this configuration.
In the last experiment (Table 2 For three challenging CBCT scans, Figure 6 presents the ground-truth and the predicted landmarks with respect to the 5-landmarks configuration DB architecture, annotated in blue and pink, respectively. We evaluated 5-landmarks configuration for both MLP and the DB architectures using variational-dropout as regularizer (Table 3). For 4-folds, we observed that DB architecture was robust and fast-to-converge. Although, the performances were statistically similar for the mandibular landmarks, this was not the case for the maxillary and the nasal bone landmarks. The performance of the MLP architecture degrades notably compared to the decrease in the DB architecture for the maxilla and nasal bone landmarks. (Table 3)  In comparison of 5-landmarks and 6-landmarks configurations (Table 3), we observed that 5-landmarks configuration is good at capturing the relations on the same bone. In contrast, 6landmarks configuration was good at capturing the relations on the neighbouring bones. Although, the error rates were less than 2mm, potentially redundant information induced by the N a landmark in the 6-landmarks configuration caused the performance to decrease notably for the mandibular landmarks compared to the 5-landmarks configuration.

3-landmarks and 5-landmarks configurations
9-landmarks configuration performed statistically similar to 5-landmarks configuration, however, due to 81 RUs employed for the 9-landmarks, the training was slower.
Although direct comparison was not possible, we compared our results with Gupta et al. 10 based on the landmark distances. We found that our results were significantly better for all landmarks except the N a landmark. The framework proposed at 10 uses an initial seed point using a 3D template registration at the inferior-anterior region where fractures are the most common. Eventually, any anatomical deformity that alters the anterior mandible may cause an error in the seed localization which can lead to a sub-optimal outcome.
We evaluated the performance of the proposed system when variational 26 and targeted 27 dropouts were employed. Although there was no accuracy-wise statistically significant difference between dropouts, convergence of the systems were relatively fast (around 20 epochs compared to 100 when using regular dropout) for the MLP architecture. Hence, for the MLP architecture, in terms of computational resources, variational and targeted dropout implementations were far more efficient for our proposed system. This is particularly important because when there are a large number of RUs, one may focus more on the efficiency rather than accuracy. When the DB architecture was employed, we did not observe any performance improvement among different dropout implementations.

Discussions and Conclusion
We proposed the RRN framework which learns spatial dependencies between CMF landmarks in an end-to-end manner. , and again no deformation or pathology presence exist therein. In contrast to these methods, our method considers a very small area as landmark, and we use extremely challenging pathological cases, which also differentiates the current work from our previous work where we used a segmentation-based approach in the geodesic space.
Our relational reasoning framework, which is a model-based approach, can generalize well to the unseen data. Hence, once trained, RRN can be used at the same testing precision to detect the missing landmarks of the unseen data taken at completely different conditions. This would afford better outcomes for precision medicine and complex CMF deformities. In our experiments, we first evaluated this claim using a dataset with a high amount of bone deformities in addition to other CBCT challenges. We observed that (1) despite the large amount of deformities that may exist in the CMF anatomy, there is a functional relation between the CMF landmarks, and (2) RNN frameworks are strong enough to reveal this latent relation information. Next, we evaluated the detection performance of five different configurations of the input landmarks to find out the optimum configuration. We observed that not all landmarks are equally informative in the detection performance. Some landmark configurations are good in capturing the local information, while some have both good local and global prediction performance. Overall, per-landmark error for the 6-landmarks configuration is less than 2mm, which is considered as a clinically acceptable level of success.
In our implementation, we showed that other deep-learning networks can be integrated well into our platform as long as features are encoded via RUs. While one may argue whether changing specific parameters could make these predictions better. However, such incremental explorations are kept outside the main paper but worth exploring in future studies from an optimization point of view. Moreover, for now RRN only employ spatial information (proof-of-concept stage), its extension could include using shape space learned landmark relationships as a conditional shape prior. Similarly, the use of those learned relationship as a look up table (atlas) is another venue that needs further exploration.
Our study has a number of limitations. For instance, we confined ourselves to manifold data (position of the landmarks and their geometric relations) without use of appearance information because one of our aims was to avoid explicit segmentation to be able to use simple geometric reasoning networks. As an extension of this study, we will incorporate appearance features from medical images to explore whether these features are superior to purely geometric features, or combined (hybrid) features can have additive value in this research domain. One alternative way to pursue the research that we initiated herein will be to explore deeper and more efficient networks.
Hence explore how to scale up in to a much wider platform where large number of landmarks and various clinical problems are addressed. We believe that such advances will improve the current technology for 3D visualization and even afford embedding augmented reality to treatment and surgical planning.

Disclosures
No Conflict of Interest.

Acknowledgments
We thank Mary McIntosh for helping data collection and landmarking. This work is partially supported by the NIH grants: R01-CA246704 and R01-CA240639.

Data, Materials, and Code Availability
Data was collected under IRB approved (PI: Janice S. Lee), and maybe accessible under certain agreement with the NIDCR. Code is available upon request.