Automatic Registration between Cone-Beam CT and Scanned Surface via Deep-Pose Regression Neural Networks and Clustered Similarities

TL;DR

This paper introduces an automatic registration method for cone-beam CT and scanned dental models using deep neural networks and clustered similarities, achieving high accuracy without human intervention or prior information.

Contribution

The study presents a novel fully automatic registration approach combining deep-pose regression neural networks with cluster-based optimization, improving accuracy and robustness in dental imaging alignment.

Findings

Registration accuracy improved by 30.77% to 70% over existing methods.

Method does not require human interaction or prior information.

Effective removal of outliers in iso-surface regions.

Abstract

Computerized registration between maxillofacial cone-beam computed tomography (CT) images and a scanned dental model is an essential prerequisite in surgical planning for dental implants or orthognathic surgery. We propose a novel method that performs fully automatic registration between a cone-beam CT image and an optically scanned model. To build a robust and automatic initial registration method, our method applies deep-pose regression neural networks in a reduced domain (i.e., 2-dimensional image). Subsequently, fine registration is performed via optimal clusters. Majority voting system achieves globally optimal transformations while each cluster attempts to optimize local transformation parameters. The coherency of clusters determines their candidacy for the optimal cluster set. The outlying regions in the iso-surface are effectively removed based on the consensus among the optimal clusters. The accuracy of registration was evaluated by the Euclidean distance of 10 landmarks on a scanned model which were annotated by the experts in the field. The experiments show that the proposed method's registration accuracy, measured in landmark distance, outperforms other existing methods by 30.77% to 70%. In addition to achieving high accuracy, our proposed method requires neither human-interactions nor priors (e.g., iso-surface extraction). The main significance of our study is twofold: 1) the employment of light-weighted neural networks which indicates the applicability of neural network in extracting pose cues that can be easily obtained and 2) the introduction of an optimal cluster-based registration method that can avoid metal artifacts during the matching procedures.

Figures

Tables

Table 1. TABLE I: MEAN AND STANDARD DEVIATION OF EUCLIDEAN DISTANCE ERROR FOR 8-FOLD CROSS-VALIDATION OF THE POINT AND LINE REGRESSIONS

Input	Point [mm]	Line [°]
$I_p$	$2.08\pm 1.45$	$2.54\pm 2.56$
$M_p$	$8.01\pm 6.20$	$2.59\pm 2.37$

Table 2. TABLE II: THE MEAN EUCLIDEAN DISTANCE ERRORS OF LANDMARKS AND RUN-TIME

Methods	Mean [mm]	Duration [s]
Three-points	$2.73\pm 0.51$	N/A
Three-points + ICP [16]	$6.29\pm 2.71$	$1.37\pm 0.40$
Three-points + Sparse ICP [18]	$3.11\pm 0.44$	$23.11\pm 4.14$
Three-points + Go-ICP [49]	$6.30\pm 2.68$	$20.74\pm 3.32$
Ours (only initial)	$4.39\pm 1.48$	$3.67\pm 0.06$
Ours	$\text{1.89}\pm \text{0.82}$	$17.68\pm 7.76$
Ours (w/o stochastic)	$\text{2.37}\pm \text{0.56}$	$7.14\pm 1.99$

Equations

$$C{v}_i={\lambda }_i{v}_i,i\in \{0,1,2\},$$

$$\chi (X,\mathbf{p},\theta ;W)={||\mathbf{p}-{\mathbf{y}}_1||}_2+\alpha {||\theta -y_2||}_2+\beta {||W||}_2^2$$

$$L(M,I;C^i,T_{{\theta }_i})=-\sum _{\mathbf{v}\in C^i}<T_{{\theta }_i}^{-1T}(\overrightarrow{M(\text{v})}),\nabla I(T_{{\theta }_i}(\text{v}))>,$$

$$\underset{{\theta }_i}{\mathrm{argmin}}-L(M,I;C^i,T_{{\theta }_i}).$$

$$d(C;T_1,T_2)=\frac{1}{|C|}\sum _{\mathbf{v}\in C}{||T_1(\mathbf{v})-T_2(\mathbf{v})||}_2,$$

$$e(C^i,C^j)=d(C^j;T_i,T_j).$$

$$c_i=e(C^i,C_0^{\ast })+e(C^i,C_1^{\ast }),$$

$$C^{min}=\underset{C^i}{\mathrm{argmax}}{c}_i.$$

$${\theta }^{\ast }=\underset{\theta }{\mathrm{argmin}}-L(M,I;C^{\ast },T_{\theta }).$$

$$d({\mathbf{l}}_i;T_{{\theta }^{\ast }})=\frac{1}{10}\sum _{i=0}^{10}{||{\mathbf{l}}_i-T_{{\theta }^{\ast }}({\mathbf{l}}_i)||}_2,$$

Full text

I Introduction

Computerized registration of medical images is a challenging problem for many clinical applications (e.g., follow-up studies, surgical planning, and augmented reality [1, 2, 3, 4]). In dental application, registration between maxillofacial cone-beam computed tomography (CT) images and a scanned model is an essential prerequisite in surgical planning for dental implants or orthognathic surgery [5, 6, 7, 8, 9, 10, 11, 12, 6, 3, 13]. Rigid transformation, which has a relatively small number of parameters, can be a viable method in the application. However, the registration problem between two different imaging protocols (i.e., optically scanned surface model and CT image) is a challenging task [6]. Moreover, the metal artifacts in cone-beam CT images hinder the accuracy of the registration [14]. Manual registration or initial registration between CT and scanned dental model is a time-consuming task. Most of the clinical applications require the three points for the initial registration. The term “three-points” refers to the registration using manual three points depiction. Thus, even the state-of-the-art clinical applications relying on the three-points step cannot be fully automated.

Many methods have been proposed to address the rigid registration problem in medical imaging applications. Intensity-based methods attempt to attain optimal registration parameters (i.e., transformation matrix) by intensity similarity criteria [15]. On the other hand, feature-based methods attempt to extract useful features (e.g., spatial points or descriptors) to match correspondences [15]. Iterative closest points (ICP) method [16, 17, 18] is widely used to perform point-to-point registrations. In dental field, the ICP method is an effective approach that can be applied to many applications such as real-time registration of optical scanning [19]. However, owing to the different protocols of the target and moving images, the ICP method is not directly applicable to the CT-model registration task. Some studies have been proposed to extract corresponding surface points from CT to apply ICP method with and without features [20, 21, 22, 23, 24, 6, 25]. Markers (i.e., manual landmarks) have also been used to perform registration between CT and optically scanned dental model [26, 27, 28, 29, 30, 31, 14]. It is also clearly known that registration performance can be increased by confident region priors [32]. The major limitations of the previously proposed methods are twofold: 1) initial registration is either not automated or lacks robustness and 2) the performance of matching is seriously affected by the extent of metal artifacts. The performance cannot guarantee robustness in the majority of cases due to the presence of metal artifacts.

In recent years, many studies have used deep neural network to resolve automatic registration problems in the field [33, 34, 35, 36, 37, 38, 39]. Some proposed unsupervised learning framework with certain similarity criteria (e.g., intensity and landmarks) and others proposed supervised [35, 36] or weakly supervised neural network [39] to perform registration. However, most of the previous works rely on intensity-based similarity or suffer from the deficiency of annotated data. The lack of ground-truth data makes it hard to train the neural network in a fully supervised manner. As aforementioned, although deep neural networks have been actively applied to various dental applications [40], the manual depiction of three-points is still required for registration in the clinics.

In this paper, we propose a novel and robust method to perform fully automatic registration between cone-beam CT images and the scanned dental model without any fragile priors (e.g., iso-surface extractions and ICP). We first extract the alignment cues (i.e., pose) from the given CT and model via deep convolutional neural networks (CNNs) followed by the rough model alignment. Finally, fine matching is performed by optimal clusters that are obtained according to the similarity measurement. The key achievements of our method are: fully automated initial alignment and robust fine registration results independent of the presence of metal artifacts in the CT images.

The remainder of this paper is structured as follows. In Section \Romannum2, we will describe our proposed method in detail. Section \Romannum3, \Romannum4, and \Romannum5 illustrate the experimental results, discussion, and conclusion, respectively.

II Methodology

Our method consists of two steps: 1) deep-pose regression and 2) optimal cluster-based matching (Fig. 1). The first step aligns the model to CT in rough conditions (i.e., initial registration). The second matching procedure performs fine registration via optimal clusters. The details of each method are described in the following sections.

II-A Deep-Pose Regression

In this step, initial registration between CT and the model is performed. Directly using 3-dimensional (3D) information is too complex and noise-dependent. Therefore, we reduced the dimension to 2D for a robust initial alignment. Let CT image $I(\textbf{x})\longrightarrow\mathbb{R}$, where $\textbf{x}\in\Omega$ ($\Omega\in\mathbb{R}^{3}$). Maximum intensity projection (MIP) image with respect to x-axis direction, $I_{p}$, is generated from $I$ (indicated by yellow arrows in Fig. 1a). As for the scanned model, a synthetic depth image is generated for the main axis (i.e., full-arch visible axis). The main axis can be easily obtained by principal component analysis (PCA) [41, 42]. In a given triangular mesh model $M=\{V,E\}$, where $V$ and $E$ are sets of vertices and edges, let a 3D vector $\textbf{v}\in{V}$ be a positional vector in set $V$. Defining the mean vector by $u=\frac{1}{|V|}\sum_{i}v_{i}$, the covariance matrix of an input data can be defined by $C=\frac{1}{|V|}\sum_{i}(v_{i}-u)(v_{i}^{T}-u^{T})$. PCA analysis can be subsequently performed via eigen-decomposition or singular value decomposition [41, 43, 44]:

[TABLE]

where $\lambda_{i}$ are the eigenvalues of $C$ ($\lambda_{0}>\lambda_{1}>\lambda_{2}$). Depth image, $M_{p}$, is then generated by projecting all vertices to a tight bounding plane that has $v_{2}$ as a normal vector. Images are normalized to the range of [0..1]. Then we use the trained CNN models to acquire corresponding point and line pairs for each image (Fig. 1a). The point and line pairs are subsequently reconstructed (i.e., positioned) in the 3D domain. Since the projected bounding plane, with respect to the scanned model, is originally defined in a 3D domain, the points and lines in $M_{p}$ are automatically positioned in the 3-dimensional space. In case of $I_{p}$, the x-coordinate is set to $w/2$ (i.e., the center point for x-axis) where $w$ is the width of an input CT image. Finally, the model is spatially matched to the 3D CT image with 3D points and lines (i.e., overlapping the bounding plane of the model and a corresponding plane in a CT which is sliced by a vector in $I_{p}$ (Fig. 1a)). Whether the scanned model is the upper jaw (i.e., maxilla) or the lower jaw (i.e., mandible) is given as a prior.

For training, we used $200$ scanned models and $100$ CT images. We manually annotated (i.e., point and angle of the line) images for $M_{p}$ and $I_{p}$. The overall loss is formulated as follows.

[TABLE]

where $X$, p, and $\theta$ are the input 2D image, ground-truth 2D point and angle of the line, respectively. $W$ represents the weights of the network, $\mathbf{y_{1}}$ and $y_{2}$ are the network outputs. The network is trained according to the weighting parameters $\alpha$ and $\beta$. For training and inference, we used the traditional VGG-16 network [45] with a slight modification in the final layer to output a 3D tensor (i.e., a 2D point and an angle). Two identical neural networks were used to train each projection image (i.e., $M_{p}$ and $I_{p}$). The only difference is the final output tensor which is 6D for $I_{p}$ (i.e., a pair of a 2D point and an angle). ‘Xavier’ initialization [46] was used for initializing all the weights of the network. While training the network, we fixed the loss parameter as $\alpha=\beta=0.1$. We used the Adam optimizer [47] with batch size 64 and set the learning rate to 0.001. We decayed the learning rate by multiplying 0.1 for every 20 epochs. We trained the network for 100 epochs using an Intel i7-7700K desktop system with 4.2 GHz processor, 32 GB of memory, and Nvidia Titan Xp GPU machine. It took 1h to complete all the training procedures.

II-B Optimal Cluster-based Matching

Cluster $C_{v,r}\in{M}$ is defined by the center vertex $v$ and radius $r\in{\mathbb{R}}$, i.e., connected local vertices in $M$ constrained by $r$ centered at $v$. We use the notation $C$ in the remaining text for simplicity. Multiple clusters are automatically generated after initial registration (Fig. 1b). We limited the minimum distance among the center of clusters to $r$ (i.e., allowing overlaps). Since we already know a projection vector of a given mesh model, we prioritized the vertices that are located in the crown region (i.e., lower values in $M_{p}$ image (Fig. 3b)). Finally, we additionally added clusters in a stochastic manner. We added three clusters for each cluster by randomly positioning the center distanced less than $r$mm with $\pm 10^{\circ}$ random rotations for each axis. The presented stochastic procedure enhances the stability of the results by means of improving the probability of the accurate local matching performance.

Local optimization (i.e., registration) is performed for each $i^{th}$ cluster ($i\in[0..N-1]$, where $N$ is the number of clusters) according to the vector alignment-based similarity criteria:

[TABLE]

where $\overrightarrow{M(\textbf{v})}$ is a normal vector at vertex v, $\nabla I(\textbf{v})$ is a gradient at $\textbf{v}\in\Omega$, and $T_{\theta_{i}}$ is a transformation with respect to the parameters, $\theta_{i}$. We used downhill simplex method [48] for local cluster optimization until convergence within six dimensional parameters (i.e., three dimensional parameters for each translation and rotation; $\theta_{i}$):

[TABLE]

After performing local optimizations of the clusters, we have scattered transformation matrices for each cluster. Let $T_{i}$ represent the locally optimized global transformation matrix of a cluster, $C^{i}$. Defining the distance function of a cluster with respect to two different transformations as

[TABLE]

the mutual coherency error for two clusters’ transformation is approximated by the following equation:

[TABLE]

The final coherency error for a given cluster is defined by

[TABLE]

where $C^{*}_{0}$ and $C^{*}_{1}$ represent the first and second minimum coherency error-valued clusters with respect to $C^{i}$ (Fig. 2a). The outlying cluster is defined as

[TABLE]

The schematic view of coherency evaluation is presented in Fig. 2. Subsequently, we remove the worst coherently transformed cluster. The iteration continues until the optimal cluster set contains three clusters. The detailed procedures are illustrated in Algorithm 1. The final three optimally chosen clusters, $C^{*}$, again iterate over the fine optimization process. In the fine registration stage, we integrate the three optimal clusters and iterate over the optimization process similar to (4):

[TABLE]

III Experiments

In this section, we present a clear visualization of deep-pose regressions and initial registration results together with the optimal clusters. Accuracy and time complexity were the primary criteria for performance evaluation. Accuracy is evaluated with the ground-truth registration result obtained by the clinical experts in the field.

We acquired data from 145 subjects, each data sample including both a CT image and the scanned models of the subject’s upper/lower jaw (i.e., 290 pairs for possible registration). The images are sourced from four different multi-centers. The model is acquired by optically scanning the surface from the cast of the patients. Each optically scanned surface formed a triangular mesh structure with $35-70$ thousand triangles and $18-40$ thousand vertices. In the CT dataset, the thickness of the slices ranged from 0.2 to 2.0mm and pixel sizes ranged from 0.2 to 1.0mm. We used 100 subjects for training, and 45 subjects for testing.

III-A Deep-Pose Regressions

In this section, we visualize the point and line regressions and the corresponding initial registration results. Fig. 3 shows each result along with its corresponding $I_{p}$ and $M_{p}$. The initial registration results using points and lines are clearly shown in Fig. 3d. For quantitative analysis, we used 8-fold cross validation metric for 100 subjects. Table I shows the Euclidean distance errors for each point and line. The result shows that there is no significant variation in the inference (i.e., testing). Even if the initial deep-pose regression resulted in a slight misalignment, the following optimal cluster-based fine registration step can complement the transformation biases.

III-B Optimal Clusters

The cluster-based fine registration algorithm is visualized in Fig. 4. As shown in Algorithm 1, crown-prioritized clusters are generated for local optimization, and the three optimal clusters that have the most coherent local transformations are selected (Fig. 4c). The final registration procedure harnesses the three remaining optimal clusters as demonstrated in (9). The metal artifact regions whose crown boundaries are difficult to delineate are automatically neglected by the non-coherent transformation (Fig. 4). Consequently, the clusters that have relatively clear boundaries in the corresponding CT image were used for the final registration. We used $r=10$mm for all the experiments. The average number of initial clusters was 17, resulting in 68 clusters in total, including the stochastic clusters.

III-C Performance Evaluation

The evaluation results were compared with the ground-truth registration results provided by the experts in the field. In assessing the accuracy of the results, we used the distance between landmarks on surface models, which were also marked by the experts. Table II shows the mean Euclidean distance errors of landmarks and the elapsed time for performing registration on 45 test subjects (i.e., 90 registration results in total). Distance errors were calculated based on the ten landmark points:

[TABLE]

where $\mathbf{l}_{i}$ is a landmark point on the surface model and $T_{\theta^{*}}$ is the optimal transformation. We compared our proposed algorithm with three-points manual registration followed by ICP-based methods [16, 17, 49, 50, 18]. The methods are state-of-the-art approaches that are frequently used in industry and many types of research. To obtain the target points from a CT image, we first extracted iso-surface of the crown with the threshold value $1000$. It is well known that the bone structure can be segmented by $1000$ Hounsfield unit in CT images [51]. The three-points method is used as a base metric, and the basic ICP [16], sparse ICP [18], and Go-ICP [49] methods are used in the subsequent fine registration. We also experimented further by eliminating the stochastic clusters from our method.

The result indicates that our proposed method achieved the highest accuracy among all the state-of-the-art methods. The ICP-based methods showed poor accuracy compared to the three-points manual registration that was carefully performed by the experts. The degradation of the accuracy from ICP-based methods is reasonable because the iso-surfaces extracted from CT images were severely affected by noise (Fig. 10). Our method, on the other hand, successfully achieved the highest accuracy compared to other methods. The aforementioned experimental variant that skips the additional generation of clusters (i.e., proceeds without stochastic addition of clusters) also outperformed the ICP-based methods. Since the randomly added clusters have a higher chance of avoiding the metal artifact regions, the stochastic method improved the registration accuracy of our method. Stochastic rotations also aid the clusters to better match with the model in the initial state. The analysis of the local regions (i.e., local optimization and coherency analysis) makes our method relatively slower than the basic ICP method. However, unlike other methods, the complete automatization of the procedure and superior accuracy may compensate for the time complexity in dental applications.

Some of the samples are visualized in Fig. 10. Fig. 10a shows manual registration with three-points. The first row in Fig. 10 clearly demonstrates that all the methods accurately perform matching when there is no metal artifact in the CT image. In the case of the crown regions with metal artifact, the ICP-based methods fail to precisely match the regions even when they are provided with the three-points as a prior (the second and third rows in Fig. 10b$\sim$d). The critical drawback of the ICP-based methods is that noisy region (i.e., metal artifact) strongly affects the target (i.e., iso-surface from CT image), and there is no explicit method to remove the metal noise in the iso-surface. Our proposed method (Fig. 10e), on the other hand, robustly matched the scanned models to CT images in every condition. The main factor of superior accuracy is the selective use of regions unaffected by noise and exclusion of noisy regions from the scope of consideration in the subsequent process (Fig. 4).

III-D Parameter Study

In this section, we study the radius parameter $r$ of a cluster. The size of the cluster affects the accuracy, run-time, and the number of clusters. Fig. 11 shows the relationship between the radius, accuracy, and run-time. The optimum radius, which attained the highest accuracy, was 10mm (Fig. 11). The accuracy decreased steeply with bigger radius. It demonstrates that bigger clusters have a higher probability of containing noisy regions. The noise may either be metal artifacts or irrelevant portions that do not contain any corresponding features to the ones that exist in a CT scan. The duration time monotonically decreases as the radius increases from the optimum while losing the accuracy. It is obvious that a bigger radius leads to a decrease in the number of clusters and vice versa. The most time-consuming procedure of our proposed algorithm is the optimal cluster-based matching (i.e., local optimization). Local optimization with many small clusters is computationally costly (Fig. 11). In addition, iterative coherency analysis of transformation (7) is another cause of performance bottlenecks.

The cluster-based registration can be viewed as a special case of the ICP method. We iteratively removed outlying points (i.e., cluster) and performed correspondence matching via optimal vector alignments. We implicitly applied noise removal procedures by removing the clusters. The radius parameter does affect the performance of the proposed method. The results showed stable accuracy when using 10mm radius. The stability of the algorithm is mainly achieved by taking many cluster candidates into account (i.e., coherency evaluation in (7)).

IV Discussion and Conclusion

Dental cone-beam CT image analysis is a challenging problem because of the large variety of cases and the presence of metal artifacts. To overcome the taxing manual labor of the registration between CT images and optically scanned models, we propose a fully automated CNN-based deep-pose estimations followed by cluster-based matching. Our work suggests the applicability of a neural network in simplifying the manual tasks involved in registration. Applying pose regressions with neural networks allowed the use of traditional methods in the following steps. As a result, we could avoid the pitfall of handling the entire registration task with complex neural networks. In fact, it is not feasible to solve the 3D registration problem directly with a neural network. The main limitations arise from the difficulties of generating the ground-truth data and unclear similarity criteria between images of different modality. As for fine registration, optimal cluster-based similarity matching was performed to accomplish confident region matching. Our method reduced the distance error by 30.77% to 70% when compared to other state-of-the-art registration methods. The primary factor for such improvement was the optimal local clusters that helped rule out the regions affected by metal artifacts.

The proposed algorithm is applicable to full-arch scanned models in a clinical practice. The performance of the proposed initial registration by deep regression networks shows results that are stable enough to be used for further fine registration process. Although working with full-arch scanned model is a common practice, the occasional use of partially scanned models may be problematic for the algorithm. This is because we employed the PCA method to obtain the axis for projection in the scanned model.