1.1.

Background and Related Work

1.2.

Summary of Our Contributions

The rest of this paper is organized as follows: we introduce our novel methodology and its details in Sec. 2. In Sec. 3, we present experiments and results and then we conclude the paper by discussing strengths and limitations of our study in Sec. 4.

2.1.

Overview and Preliminaries

The most frequently deformed or injured CMF bone is the lower jawbone, or mandible, which is the only mobile CMF bone.²⁸ In our previous study,⁵ 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 to the following important questions:

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:

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. Toward this goal, we propose to learn a relationship between the anatomical landmarks in two stages (illustrated in Fig. 3) based on relational units (RUs). The first stage is shown as the function $g$, which learns the pairwise (local) relations. The second stage is shown as a function $f$, which combines pairwise relations ($g$) of the first stage into a global relation based on RUs.

Fig. 3

Overview of the proposed RRN architecture: for a few given input landmarks, RRN utilizes both pairwise and combination of all pairwise relations to predict the remaining landmarks.

Figure 4 shows an example of pairwise relations for different pairs of mandible landmarks. There are five sparsely localized landmarks. The basis/reference is chosen as menton (Me), in this example, hence, four pairwise relations are illustrated from Figs. 4(a)–4(d). Figure 4(e) illustrates combined relations Figs. 4(a)–4(d) of the landmark menton (reference) with respect to other four landmarks on the mandible.

Fig. 4

For the input domain $L_{\text{input}}=\{\text{Me},{\mathrm{Cd}}_L,{\mathrm{Cor}}_L,{\mathrm{Cd}}_R,{\mathrm{Cor}}_R\}$, (a)–(d) pairwise relations of landmark menton (Me): (a) menton-condylar left, (b) menton-coronoid left, (c) menton-condylar right, (d) menton-coronoid right, and (e) combined relations of menton.

2.2.

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 clinical problem is an open problem. While anatomical similarities at the local and global levels could serve toward viable solutions, thus far, features that can represent anatomical landmarks from the medical images have not achieved the desired efficacy and interpretation.^2,29–31

We propose a new framework called relational reasoning network (RRN) to learn local and global relations of anatomical landmarks ($o_i$) through its units called RU (relationship unit). The proposed RRN architecture and its RU subarchitectures are summarized in Fig. 5. The relation between two landmarks is encoded via major spatial properties of the landmarks. We explore two architectures as RU: first one is a simple multilayer perceptron (MLP) (Fig. 5-bottom left) (similar to Ref. 16), the other one is more advanced architecture composed of dense-blocks (DBs) (Fig. 5 bottom middle). Both architectures are relatively simple compared to very dense complex deep-learning architectures. The rationale is simple when there is a less data (i.e., pairwise relation), it is natural to choose fully connected layers to keep the full spectrum of the data at hand. Similarly, when more pairwise data are available for exploration of the more complex relation, it is natural to move into convolutional operation from fully connected layers to keep dominant information while reducing the redundant information and providing a computational feasibility. Our objective is to locate all anatomical landmarks by inputting a few landmarks to RRN, which provides reasoning inferred from the learned relationships of landmarks and locate all other landmarks automatically.

Fig. 5

(a) RRN architecture for five-input landmarks $\mathrm{RRN}(L_{\text{input}})$: $L_{\text{input}}=\{\text{Me},{\mathrm{Cor}}_L,{\mathrm{Cor}}_R,{\mathrm{Cd}}_L,{\mathrm{Cd}}_R\}$, $\widehat{L}=\{\mathrm{Gn},\mathrm{Pg},\mathrm{B},\mathrm{Id},\mathrm{Ans},\mathrm{A},\mathrm{Pr},\mathrm{Pns},\mathrm{Na}\}$ and $\mu$ is the average operator. (b) Content of the pairwise reasoning block, RU of the RRN. (c) RU composed of two DBs, convolution and concatenation (C) units. (d) DB architecture composed of four layers and concatenation layers.

In the pairwise learning/reasoning stage (stage 1), five-landmarks-based system is assumed as an example network (other configurations are possible too, see experiments and results section). Sparsely spaced landmarks [Fig. 4(e)] and their pairwise relationships are learned in this stage ($g_{\theta }$). These pairwise relationship(s) are later combined in a separate DB setting in ($f_{\varphi }$). It should be 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\times (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_{\text{input}}$ and $\widehat{L}$ indicate vectors of input and output anatomical landmarks, respectively. Then, two stages of the RRN of the input domain landmarks $L_{\text{input}}$ can be defined as

2.3.

Pairwise Relation ($g_{\theta }$)

For a given a few input landmarks ($L_{\text{input}}$), our objective is to predict the 3D spatial locations of the target domain landmarks ($\in \widehat{L}$) by using the 3D spatial locations of the input domain landmarks ($\in L_{\text{input}}$). With respect to relative locations of the input domain landmarks, we reason about the locations of the target domain landmarks. The RU function $g_{\theta }(o_i,o_j)$ represents the relation of two input domain landmarks $o_i$ and $o_j$ where $i\ne j$ [Figs. 4(a)–4(d)]. The output of $g_{\theta }(o_i,o_j)$ describes relative spatial context of two landmarks, defined for each pair of input domain landmarks (pairwise relation at Fig. 5). According to each input domain landmark $o_i$, the structure of the manifold is captured through mean of all pairwise relations [represented as $G_{{\theta }_i}$ at Eq. (1)].

2.4.

Global Relation ($f_{\varphi }$)

The mean pairwise relation $G_{\theta i}$ is calculated with respect to each input domain landmark $o_i$, and it is given as input to the second stage where global (combined) relation $f_{\varphi i}$ is learned. $f_{\varphi i}$ is an RU function and the output of $f_{\varphi i}$ is the predicted 3D coordinates of the target domain landmarks ($\in \widehat{L}$). In other words, each input domain landmark $o_i$ learns and predicts the target domain landmarks by the RU function $f_{\varphi i}$. The terminal prediction of the target domain landmarks is the average of individual predictions of each input domain landmark, represented by $\mathrm{RRN}(L_{\text{input}};\theta ,\varphi )$ at Eq. (1). There are totally $n^2$ RUs in the architecture. The number of trainable parameters used for each experimental configuration are directly proportional with $n^2$ (Fig. 6). Since all pairwise relations are leveraged under $G_{{\theta }_i}$ and $f_{\varphi }$ with averaging operation, we can conclude that RRN is invariant to the order of input landmarks (i.e., permutation-invariant).

Fig. 6

Five experimental landmark configurations for experimental explorations. $L_{\text{input}}$: input landmarks and $\widehat{L}$: output landmarks, and #RUs indicate the number of relational units. Landmarks are visualized using reference standard bones for illustrative purposes; in our implementation there is no explicit segmentation exist.

2.5.

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:

2.6.

Variational Dropout

Dropout is an important regularizer employed to prevent overfitting at a cost of 2 to 3 times (on average) increase in training time.³² For efficiency reasons, speeding up dropout is critical and it can be achieved by a variational Bayesian inference on the model parameters.²⁶ 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 $\omega$ such that $y=F_{\omega }(x)$. In the Bayesian approach, given the input and output datasets $X,Y$, we seek for the posterior distribution $p(\omega |X,Y)$, by which we can predict output $y^{*}$ for a new input point $x^{*}$ by solving the integral³³

In practice, this computation involves intractable integrals.²⁶ To obtain the posterior distributions, a Gaussian prior distribution $N(0,I)$ is placed over the network weights³³ which leads to a much faster convergence.²⁶

2.7.

Targeted Dropout

Alternatively, we also propose to use targeted dropout for better convergence.²⁷ Given a neural network parameterized by $\mathrm{\Theta }$, the goal is to find the optimal parameters $W_{\mathrm{\Theta }}(.)$ such that the loss $\text{Loss}(W_{\mathrm{\Theta }})$ is minimized. For efficiency and generalization reasons, $|W_{\mathrm{\Theta }}|\le k$, only $k$ weights of highest magnitude in the network are employed. In this regard, deterministic approach is to drop the lowest $|W_{\mathrm{\Theta }}|-k$ weights. In targeted dropout, using a target rate $\gamma$ and a drop out rate $\alpha$, first a target set is generated with the lowest weights with the target rate $\gamma$. Next, weights are dropped out in a stochastic manner from the target set at a certain dropout rate $\alpha$.

2.8.

Landmark Features

Pairwise relations are learned through RU functions. Each RU accepts input features to be modeled as a pairwise relation. It is desirable to have such features characterizing landmarks and interactions with other landmarks. These input features can either be learned throughout a more complicated network design, or through feature engineering. In this study, for simplicity, we define a set of simple yet explainable geometric features. Since RUs model relations between two landmarks ($o_A$ and $o_B$), we use 3D coordinates of these landmarks (both in pixel and spherical space), their relative positions with respect to a well-defined landmark point (reference), and approximate size of the mandible. The mandible size is estimated as the distance between the maximum and the minimum coordinates of the input domain mandibular landmarks (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.

Table 1

Input landmarks have the following feature(s) to be used only in stage I. 19D feature vector includes only structural information.

3.1.

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 CBCT system (Hitachi Medical Corporation, Tokyo, Japan) was used to scan the data at 10 mA and 100 Kvp. The radiation dosage for each scan was around 300 mSv. To handle the computational cost, each patient’s scan was resampled from $512\times 512\times 512$ to $256\times 256\times 512$. In-plane resolution of the scans was noted (in mm) either as $0.754\times 0.754\times 0.377$ or $0.584\times 0.584\times 0.292$. 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.⁵ Briefly, 3% of the whole data set was including CBCT scans with extreme deformation and artifacts, while 11% of the data set was including cases with large-scale tissue or bone deformations, artifacts, or missing bones. 16% of the data set was including minor tissue deformation and/or metal or other artifacts. The remaining 70% of the data was either no visible artifacts or minor problems in visual assessment. These statistics were obtained by participating two experts, blindly to each other, and qualitatively they were asked to evaluate the scans visually.

The data was annotated independently by three expert interpreters, one from the NIH team, and two from UCF team. Interobserver agreement values were computed as $\sim 3\text{pixels}$. Experts used freely available 3D Slicer software for the annotations.⁵

3.2.

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 was 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.³⁰ 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 $\pm 5\text{pixels}$, defined empirically based on the resolution of the images as well as the observed high deformity of the bones. We generated a dataset with $\sim 100\mathrm{K}$ landmarks, which is empirically evaluated as a sufficiently large dataset.

3.3.

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.

3.4.

Input Landmark Configurations

In our experiments, there were three groups of landmarks (See Fig. 2) defined based on the bones they reside: Mandibular $L_1=\{o_1,\dots ,o_9\}$, Maxillary $L_2=\{o_{10},\dots ,o_{13}\}$, and Nasal $L_3=\{o_{14}\}$, where subscripts in $o$ denote the specific landmark in that bone:

In each experiment, as detailed in Fig. 6, we designed a specific input set $L_{\text{input}}$ where $L_{\text{input}}\subseteq L_1\cup L_2$, $|L_{\mathrm{input}}|=n$ and $1<n\&lt;=(|L_1|+|L_2|)$. The target domain landmarks for each experiment were $\widehat{L}=(L_1\cup L_2\cup L_3)\setminus L_{\text{input}}$ and $|\widehat{L}|=m$ such that $n+m=14$. With carefully designed input domain configurations $L_{\text{input}}$, and pairwise relationships of the landmarks in the input set, we seek the answers to the following questions previously defined as Q1 and Q2 in Sec. 2:

Overall, we designed five different input landmark configurations called three-landmarks regular, three-landmarks cross, five-landmarks, six-landmarks, and nine-landmarks (Fig. 6). Each configuration is explained in Sec. 3.6.

3.5.

Training

The MLP RU was composed of three fully connected layers, two batch normalizations, and two ReLUs (Fig. 5). The DB RU architecture contained 1 DB, which was composed of four layers with a growth rate of 4. We used a batch size of 64 for all experiments. For the five-landmarks configuration, there were 6,596,745 and 11,068,655 trainable parameters for the MLP and the DB architectures, respectively. We trained the network for 100 epochs on 1 NVIDIA Titan-XP GPU with 12 GB memory using the MLP architecture with the regular dropout compared to 20 epochs with the variational and targeted dropout implementations. For the DB architecture, it converged in around 20 epochs independent of the dropout implementation employed.

3.6.

Experiments and Results

We ran a set of experiments to evaluate the performance of the proposed system using a fourfold cross-validation. We summarized the experimental configurations in Fig. 6, error rates in Table 2, and corresponding renderings in Fig. 7. The method achieving the minimum error for a corresponding landmark is colored the same as the corresponding landmark at Table 2. As shown by results, the minimum number of the input landmarks required for successful identification of other landmarks is determined as 3.

Fig. 7

Landmark annotations using the five-landmarks configuration: Ground truth in blue and computed landmarks in pink. (a) Genioplasty/chin advancement (male 43 year old), (b) malocclusion (mandibular hyperplasia, maxillary hypoplasia) surgery (male 19 year old), (c) malocclusion (mandibular hyperplasia, maxillary hypoplasia) surgery (female 14 year old). Note that landmarks are shown on the volume-rendered CBCT scans; there is no segmentation conducted.

Table 2

Landmark localization errors (mm). The symbol “—” means not applicable (N/A). Bold values represent the best performance obtained from the ablation and benchmarking experiments.

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 five-landmark experiment (See Table 2).

In the first experiment (Fig. 6, 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.⁵ We named our first configuration as “five-landmarks” where closely spaced maxillary and nasal bones’ landmarks are predicted based on the relation of sparsely spaced landmarks (Fig. 6). In the five-landmarks RRN architecture, there were totally 25 RUs. In the second experiment (Fig. 6, Experiment 2), we explored the impact of a configuration with a smaller number of input mandibular landmarks on the learning performance. Compared to the five sparsely spaced input landmarks, we learned the relation of three landmarks, Me, ${\mathrm{Cd}}_L$, and ${\mathrm{Cd}}_R$, and predicted the closely-spaced landmark locations (as in the five-landmarks experiment) plus superior-anterior landmarks ${\mathrm{Cor}}_L$ and ${\mathrm{Cor}}_R$ and maxillary and nasal bones’ landmark locations. The network was composed of nine RUs. The training was relatively fast compared to the five-landmarks configuration due to small number of RUs. We named this method as “three-landmarks regular.”

After observing statistically similar accuracy compared to the five-landmarks method for the closely-spaced landmarks ($P>0.005$), and high error rates at the superior–anterior landmarks ${\mathrm{Cor}}_L$ and ${\mathrm{Cor}}_R$, we set up a new experiment, which we named “three-landmarks cross” (Fig. 6, Experiment 3). In this configuration, the third experiment, we used one superior–posterior and one superior–anterior landmarks on the right and left sides, respectively. This network was similar to three-landmarks regular one in terms of number of RUs used.

In the fourth experiment (Fig. 6, Experiment 4), we evaluated the performance of the system in learning the closely-spaced mandibular landmarks (Gn, Pg, B, Id) and the maxillary landmarks (ANS, A, Pr, PNS) using the relation information of the sparsely spaced and the nasal-bones landmarks which is named as “six-landmarks.” There are a total of 36 RUs in this configuration.

In the last experiment (Fig. 6, Experiment 5), we aimed to learn the maxillary landmarks (ANS, A, Pr, PNS) and nasal bones landmark (Na) using the relation of the mandibular networks; hence, this network configuration is called “nine-landmarks.” The architecture was composed of 81 RUs. Owing to the high number of RUs in the architecture, the training of this network was the slowest among all the experiments performed.

For three challenging CBCT scans, Fig. 7 presents the ground-truth and the predicted landmarks with respect to the five-landmarks configuration DB architecture, annotated in blue and pink, respectively. We evaluated five-landmarks configuration for both MLP and the DB architectures using variational-dropout as regularizer (Table 2). For fourfolds, 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.

Three-landmarks and five-landmarks configurations (Table 2) performed statistically similar for the mandibular landmarks. Interestingly, both three-landmarks configurations performed slightly better for the neighboring bone landmarks. This reveals the importance of optimum number of landmarks in the configuration.

In comparison of five-landmarks and six-landmarks configurations (Table 2), we observed that five-landmarks configuration is good at capturing the relations on the same bone. In contrast, six-landmarks configuration was good at capturing the relations on the neighboring bones. Although, the error rates were $<2\mathrm{mm}$, potentially redundant information induced by the Na landmark in the six-landmarks configuration caused the performance to decrease notably for the mandibular landmarks compared to the five-landmarks configuration.

Nine-landmarks configuration performed statistically similar to five-landmarks configuration, however, due to 81 RUs employed for the nine-landmarks, the training was slower.

Although direct comparison was not possible, we compared our results with Gupta et al.¹⁰ based on the landmark distances. We found that our results were significantly better for all landmarks except the Na landmark. The framework proposed at Ref. 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 suboptimal outcome.

We evaluated the performance of the proposed system when variational²⁶ and targeted²⁷ dropouts were employed. Although there was no statistically significant difference between dropouts in terms of accuracy, convergence of the systems was 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.

In landmarking, extreme performance would be very important. For example, the outliers would hamper the entire planning. Therefore, we have carefully checked the outliers for each landmarking problem, and found that there are $<10$ outliers in a total of 250 patients’ scans for each configuration, and the highest number of outliers was 7. In the best working experimental setup (six-landmark configuration), for both menton and condylar left, the highest errors we obtained with outliers were 1.5 mm. For coronoid left and right, errors of 2.75 and 0.3 mm were obtained, respectively. For infradentale, we obtained 1.5 and 1.7 mm errors by two outliers, all measured in volumetric spaces. Our results were consistent and robust to outliers, as hypothesized before. In the three-landmark configuration, the highest error we obtained with an outlier was 5 mm to detect coronoid right.

It is also important to explore whether severity of the cases or metal artifacts can influence the final results, considered under the robustness measure. We specifically evaluated the performance of our method within the subsets of data containing significant artifacts and others (30% versus 70%, see data subsection). Statistically, there was no difference in the performance of our landmarking method when applied a $t$-test between these two groups, indicating a robustness of our method.