Three-Dimensional Tooth Model Reconstruction Using Statistical Randomization-Based Particle Swarm Optimization

Abstract

The registration between images is a crucial part of the 3-D tooth reconstruction model. In this paper, we introduce a registration method using our proposed statistical randomization-based particle swarm optimization (SR-PSO) algorithm with the iterative closet point (ICP) method to find the optimal affine transform between images. The hierarchical registration is also utilized in this paper since there are several consecutive images involving in the registration. We implemented this algorithm in the scanned commercial regular-tooth and orthodontic-tooth models. The results demonstrated that the final 3-D images provided good visualization to human eyes with the mean-squared error of 7.37 micrometer² and 7.41 micrometer² for both models, respectively. From the results compared with the particle swarm optimization (PSO) algorithm with the ICP method, it can be seen that the results from the proposed algorithm are much better than those from the PSO algorithm with the ICP method.

1. Introduction

A 3-D tooth model reconstruction is a very crucial part of orthodontics. It is a helpful tool in diagnostic or examination and also problem identification in the treatment planning process in other dental caries in adults and children [1]. This is especially true in children in which dental caries are one of the most children chronic diseases [2]. Nowadays, the most powerful tools are image reconstructions from sophisticated devices, for example, CT or laser. There are several research works on CT image reconstruction [3] and the reconstruction from multimodal images [4,5,6]. Unfortunately, oral healthcare is not sufficient and dental care access is limited, particularly in rural areas [7]. In Thailand, the dental innovation foundation under royal patronage has provided dental care access in rural communities for a long time. However, one of the difficulties in providing dental care to children in those areas is how to take a 3-D tooth image inside their mouths. It is not easy to use any sophisticated imaging devices when the units are applied in the countryside due to the distance or the location limitation.

Therefore, a 3-D reconstruction system from optical images is needed. Still, one of the essential processes in 3-D image reconstruction is image registration. In the literature, there exist some 2-D medical image registration research works [8,9,10,11,12]. However, in our case, the 3-D registration is more suitable. There exist some 3-D medical image registration methods [13,14] that utilize several features in the registration process including point-cloud coordinates representing the 3-D shapes of objects. These coordinates have also been used in the registration process shown in [15,16,17,18,19,20]. All the mentioned research works utilize a variation of the particle swarm optimization (PSO) to find the matching location between the source and target images. Most of the existing registration methods find the matching locations/points based on rotation and/or translation only. However, in the 3-D affine transform, there are other types of transform, e.g., scaling, shearing, and reflection.

In this paper, we develop a system that can create a 3-D image from optical images. However, it is not easy to use real images taken from children due to a research ethical approval requirement. Therefore, we postulate scanned images from two commercial tooth models and then create point-cloud images. We propose the statistical randomization-based particle swarm optimization (SR-PSO) algorithm to find an appropriate affine transform (shearing, rotation, scaling and translation) between images for the registration purpose. The proposed SR-PSO algorithm is based on a modification of the particle swarm optimization [21] to cope with the premature convergence [22] and to improve exploration and exploitation of the algorithm [23]. After that, the iterative closet point (ICP) method [24,25] is used to refine the resulting registration because the ICP method has been proved in several research works [8,17] that it can help to refine registered results. Finally, the 3-D tooth models are reconstructed.

2. Registration Method

In a 3-D registration problem, the geometry transform can be calculated from the relation between two point-cloud images. The more they are related, the higher quality of the transformation. The highest quality transformation corresponds to the minimization of the following statement:

$$H^{*}=\underset{}{\arg \min }O(H(Q),P)$$

(1)

where P is the target point-cloud matrix ([p_i]_M×4, M is the number of target point-cloud points), Q is the source point-cloud matrix ([q_j]_N×4, N is the number of source point-cloud points), and H is the geometry transform. Finally, O(⋅) is an objective function. Since the transformation H is estimated by finding the nearest neighbor [26] between a set of point-pairs (p_j, q_j), the minimum error of the distance between two corresponding points can be considered [27]. Using the mean squared error (MSE), hence, the minimization problem, in this case, is calculated as

$$H^{*}=\arg {\min }_H\frac{1}{N}{{\displaystyle \sum _{j=1}^N\left({\mathit{q}}_j\cdot H^T-{\mathit{p}}_j\right)}}^2$$

(2)

$${\mathit{p}}_j=\underset{{\mathit{p}}_i\in P}{\arg \min }\Vert {\mathit{q}}_j.H^T-{\mathit{p}}_i\Vert$$

(3)

The transformation H is the 3-D transformation with 15 unknown parameters, i.e., 3 parameters from scaling (S), 3 parameters from translation (T), 3 parameters from rotation (R), and 6 parameters from shearing (SH) [28]. The matrix H is computed as

H= T × S × R × SH,

(4)

where $T=\left[\begin{array}{cccc}1& 0& 0& t_x\\ 0& 1& 0& t_y\\ 0& 0& 1& t_z\\ 0& 0& 0& 1\end{array}\right]$, $S=\left[\begin{array}{cccc}{s}_x& 0& 0& 0\\ 0& s_y& 0& 0\\ 0& 0& s_z& 0\\ 0& 0& 0& 1\end{array}\right]$, and $SH=\left[\begin{array}{cccc}1& s{h}_1& s{h}_2& 0\\ s{h}_3& 1& s{h}_4& 0\\ s{h}_5& s{h}_6& 1& 0\\ 0& 0& 0& 1\end{array}\right]$, and

$$R=\left[\begin{array}{cccc}\cos \left({\theta }_y\right)\cos \left({\theta }_z\right)& -\cos \left({\theta }_y\right)\sin \left({\theta }_z\right)& \sin \left({\theta }_y\right)& 0\\ \sin \left({\theta }_x\right)\sin \left({\theta }_y\right)\cos \left({\theta }_z\right)+\cos \left({\theta }_x\right)\sin \left({\theta }_z\right)& -\sin \left({\theta }_x\right)\sin \left({\theta }_y\right)\sin \left({\theta }_z\right)+\cos \left({\theta }_x\right)\cos \left({\theta }_z\right)& -\sin \left({\theta }_x\right)\cos \left({\theta }_y\right)& 0\\ -\cos \left({\theta }_x\right)\sin \left({\theta }_y\right)\cos \left({\theta }_z\right)+\sin \left({\theta }_x\right)\sin \left({\theta }_z\right)& \cos \left({\theta }_x\right)\sin \left({\theta }_y\right)\sin \left({\theta }_z\right)+\sin \left({\theta }_x\right)\cos \left({\theta }_z\right)& \cos \left({\theta }_x\right)\cos \left({\theta }_y\right)& 0\\ 0& 0& 0& 1\end{array}\right]$$

(5)

Hence,

$$H=\left[\begin{array}{cccc}a& d& g& l\\ b& e& i& m\\ c& f& j& n\\ p& q& r& s\end{array}\right]$$

(6)

where

$$\begin{array}{l}a=s_x(cycz)+s{h}_3{s}_x(-cysz)+s{h}_5{s}_x(sy),\\ b=s_y(sxsycz+cxsz)+s{h}_3{s}_y(-sxsysz+cxcz)+s{h}_5{s}_y(-sxcy),\\ c=s_z(-cxsycz+sxsz)+s{h}_3{s}_z(cxsysz+sxcz)+s{h}_5{s}_z(cxcy),\\ d=s{h}_1{s}_x(cycz)+s_x(-cysz)+s{h}_6{s}_x(sy),\\ e=s{h}_1{s}_y(sxsycz+cxsz)+s_y(-sxsysz+cxcz)+s{h}_6{s}_y(-sxcy),\\ f=s{h}_1{s}_z(-cxsycz+sxsz)+s_z(cxsysz+sxcz)+s{h}_6{s}_z(cxcy),\\ g=s{h}_2{s}_x(cycz)+s{h}_4{s}_x(-cysz)+s_x(sy),\\ i=s{h}_2{s}_y(sxsycz+cxsz)+s{h}_4{s}_y(-sxsysz+cxcz)+s_y(-sxcy),\\ j=s{h}_2{s}_z(-cxsycz+sxsz)+s{h}_4{s}_z(cxsysz+sxcz)+s_z(cxcy),\\ cx=\cos ({\theta }_x),cy=\cos ({\theta }_y),cz=\cos ({\theta }_z),\\ sx=\sin ({\theta }_x),sy=\sin ({\theta }_y),sz=\sin ({\theta }_z).\end{array}$$

(7)

It is worthwhile noting that a through j are non-rigid transformations resulting from the combination of scaling, shearing, and rotation properties. Meanwhile, l, m, and n are simply t_x, t_y, and t_z, respectively. Moreover, p, q, and r are set to 0 since they are perspective property values. Finally, s is always set to 1 because of the scaling factor.

The proposed statistical randomization-based particle swarm optimization (SR-PSO) algorithm described in the following section is used to find the optimal H. Each individual in the swarm has 15 dimensions. The search space is defined as shown in

Table 1. Parameters boundaries in optimization process.

Parameters	Lower Bound	Upper Bound
tx, ty, tz	–1.5 (cm)	1.5 (cm)
ɵx, ɵy, ɵz	–45 (deg)	45 (deg)
sx, sy, sz	0.8 (20% downscaling)	1.2 (20% upscaling)
sh1, sh2, sh3, sh4, sh5, sh6	–0.5 (cm)	0.5 (cm)

.

Statistical Randomization-Based PSO (SR-PSO) Algorithm

In this research, we modify the particle swarm optimization [21,29,30,31] following [22,23] to cope with the premature convergence, low accuracy, and to improve exploration and exploitation of the algorithm. In particular, we modify and combine the methods in [22] and [23] to exploit both of their advantages.

Let X = $\left\{{\mathit{x}}_j|j=1\dots N\right\}$ be a set of N particles in the swarm in d-dimensional feature space. ${\mathit{x}}_i^b$ and x^g are the individual best of the i^th particle and the global best of the swarm, respectively. The update equations for velocity and position of each particle are

$${\mathit{v}}_i\left(t+1\right)=w{\mathit{v}}_i\left(t\right)+c_1{r}_1({\mathit{x}}_i^b\left(t\right)-{\mathit{x}}_i\left(t\right))+c_2{r}_2({\mathit{x}}^g\left(t\right)-{\mathit{x}}_i\left(t\right))$$

(8)

$${\mathit{x}}_i\left(t+1\right)={\mathit{x}}_i\left(t\right)+{\mathit{v}}_i\left(t+1\right)$$

(9)

where r₁ and r₂ are randomly generated numbers from the uniform distribution within [0, 1]. c₁ and c₂ are the acceleration coefficients, and w is the inertia weight calculated by [29,30,31]

$$w(t+1)=w_{\max }-\frac{\left(w_{\max }\times t^{\alpha }\right)}{\left(T^{\alpha }\right)}$$

(10)

where w_max = χ, T is the number of iterations, α = 1 (normally 0 < α ≤ 2), and

$$\chi =\frac{2}{\left|2-\phi -\sqrt{{\phi }^2-4\phi }\right|},\phi >4$$

(11)

In our case, we also set the velocity limit to [v_min, v_max] where

$$v_{\max ,d}=k\times (u{b}_d-l{b}_d), \mathrm{for} 0.1\le k\le 1$$

(12)

$$v_{\min ,d}=-1\times v_{\max ,d}$$

(13)

In the experiment, we set k = 0.1.

We modify the method in [22] by introducing an extra intermediate particle by randomizing the particle’s position using Gaussian distribution in each dimension in order to increase the chance of premature convergence avoidance. The details of intermediate particles are as follows:

For the j^th dimension of intermediate particle $\left(x_j^{tmi}\right)$ with K being the number of particles in the swarm,

−

Calculated from an average of all individual best $\left(x_{ij}^b\right)$:

$$x_j^{tm1}=\frac{{\displaystyle {\sum }_{i=1}^K{x}_{ij}^b}}{K}$$

(14)

−

Calculated from a median of all individual best $\left(x_{ij}^b\right)$:

$$x_j^{tm2}=\underset{1\le i\le K}{\mathrm{median}}(x_{ij}^b)$$

(15)

−

Calculated by random generate number from Gaussian distribution with mean and standard deviation computed from individual best positions:

$$x_j^{tm3}={\sigma }_j\times Z+{\mu }_j, \mathrm{when} Z\sim N(0,1)$$

(16)

−

Calculated from the larger absolute value of the maximum and the minimum in that dimension:

$$x_j^{tm4}=\{\begin{array}{l}\underset{1\le i\le K}{\min }(x_{ij}^b) \mathrm{if} \left|\underset{1\le i\le K}{\min }(x_{ij}^b) \right|>\left|\underset{1\le i\le K}{\max }(x_{ij}^b) \right|\\ \underset{1\le i\le K}{\max }(x_{ij}^b) \mathrm{else}\end{array}$$

(17)

−

Calculated from the smaller absolute value of the maximum and the minimum in that dimension:

$$x_j^{tm5}=\{\begin{array}{l}\underset{1\le i\le K}{\min }(x_{ij}^b) \mathrm{if} \left|\underset{1\le i\le K}{\min }(x_{ij}^b) \right|<\left|\underset{1\le i\le K}{\max }(x_{ij}^b) \right|\\ \underset{1\le i\le K}{\max }(x_{ij}^b) \mathrm{else}\end{array}$$

(18)

If the l^th intermediate particle (x^tml) is the best particle among the intermediate particles and it is better than the global particle, then the global best particle will be replaced by x^tml.

To increase the chance of exploration and exploitation of x^g, we randomly select the individual best (${\mathit{x}}_k^b$) from all particles (instead of randomly selecting a particle except the best one as in [23]) using

$${\mathit{x}}_{select}^b=\mathrm{rand}({\mathbf{X}}^b) \mathrm{where} {\mathbf{X}}^b=\{{\mathit{x}}_1^b,\dots ,{\mathit{x}}_K^b\}$$

(19)

For each dimension ($x_{select,i}^b$), if a randomly generated number (∈ [0, 1]) is bigger than or equal to a certain thresholding value, we reinitialize $x_{select,i}^b$ by randomizing the number within [lb_i, ub_i], where lb_i and ub_i are the lower and upper bounds in the search space in the i^th dimension. That is

$$x_{select,i}^b=\{\begin{array}{l}\mathrm{rand}(l{b}_i,u{b}_i) \mathrm{if} \mathrm{rand}(0,1)>1-\frac{1}{d}\\ x_i^g\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{otherwise}\end{array}$$

(20)

In this case, we set the thresholding value to 1 − (1/d), where d is the number of dimensions of a particle. If ${\mathit{x}}_k^b$ is better than x^g, it is then updated using

$${\mathit{x}}^g=\{\begin{array}{l}{\mathit{x}}_{select}^b \mathrm{if} f\left({\mathit{x}}_{select}^b) < f\right({\mathit{x}}^g)\\ {\mathit{x}}^g\hspace{1em}\hspace{1em}\hspace{1em} \mathrm{otherwise}\end{array}$$

(21)

The SR-PSO algorithm is summarized as followings:

In our experiment, the optimal solution is the global best in the last population. To fine tune the registration results from the SR-PSO algorithm, after we find the optimized parameters, we utilize the iterative closest point algorithm (ICP) method as in [24,25] with the Nelder–Mead simplex method [32].

3. Experimental Results

The 3-D tooth model reconstruction system is shown in Figure 1. The SR-PSO algorithm is used to find the optimal transformation matrix H⁻¹ (transform from source point-cloud to target point-cloud). The ICP method is used to fine-tune the resultant H⁻¹. Finally, the 3-D tooth models are reconstructed based on the registered source and target point-clouds.

To test the algorithm, we ran the experiment of the generated original shape and its 3-D transformation as shown in Figure 2.

The transformation matrix used to transform the original shape described by the target point-cloud (Figure 2a) to the transformed shape described by the source point-cloud (Figure 2b) is

$$H=\left[\begin{array}{cccc}0.882050& -0.285362& -0.555884& -0.061153\\ 0.225174& 1.041540& 0.181496& 0.063487\\ 0.249299& -0.413927& 0.966936& -0.163016\\ 0& 0& 0& 1\end{array}\right]$$

(22)

Hence, the transform matrix from the source point-cloud (Figure 2b) to the target point-cloud (Figure 2a) is

$$H^{-1}=\left[\begin{array}{rrrr}\hfill 0.901895& \hfill 0.421703& \hfill 0.439339& \hfill 0.1\\ \hfill -0.143742& \hfill 0.826257& \hfill -0.237726& \hfill -0.1\\ \hfill -0.294063& \hfill 0.244980& \hfill 0.819157& \hfill 0.1\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$$

(23)

In the experiment, the parameters for the SR-PSO algorithm were set as shown in

Table 2. Statistical randomization-based particle swarm optimization (SR-PSO) algorithm’s parameters settings in the experiment.

Parameters		Values
Number of particles	K	100
Number of iterations	T	2000
Constriction coefficient	χ	0.7298
	α	0.5, 1.0, 1.5, 2.0
	φ	4.1
Personal learning coefficient	c1	1.4962
Global learning coefficient	c2	1.4962

. To demonstrate how the SR-PSO algorithm works in each step, we provide an example as follows:

Suppose an initialized particle is $H^{-1}=\left[\begin{array}{rrrr}\hfill 0.731187& \hfill -0.409449& \hfill 0.035758& \hfill 0.590480\\ \hfill 0.142237& \hfill 0.573472& \hfill 1.013846& \hfill 1.008356\\ \hfill -0.099053& \hfill -0.640562& \hfill 1.135987& \hfill 1.336424\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$, then after the velocity and position are updated, the transformation is $H^{-1}=\left[\begin{array}{rrrr}\hfill 0.791731& \hfill -0.221107& \hfill 0.249697& \hfill 0.479447\\ \hfill -0.022301& \hfill 0.854962& \hfill 0.732671& \hfill 0.571793\\ \hfill -0.010562& \hfill -0.459633& \hfill 1.176127& \hfill 0.893047\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$. The individual best of this particle is $H_b^{-1}=\left[\begin{array}{rrrr}\hfill 0.731187& \hfill -0.409449& \hfill 0.035758& \hfill 0.590480\\ \hfill 0.142237& \hfill 0.573472& \hfill 1.013846& \hfill 1.008356\\ \hfill -0.099053& \hfill -0.640562& \hfill 1.135987& \hfill 1.336424\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$ and the global best in this iteration is $H_g^{-1}=\left[\begin{array}{rrrr}\hfill 1.004627& \hfill 0.362394& \hfill -0.003666& \hfill 0.024357\\ \hfill -0.764060& \hfill 0.540426& \hfill -0.902233& \hfill 1.839488\\ \hfill 0.567330& \hfill 0.659548& \hfill 1.130398& \hfill -0.219057\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$. Then, all five intermediate particles are $H_{tm1}^{-1}=\left[\begin{array}{rrrr}\hfill 0.995780& \hfill 0.115222& \hfill -0.006401& \hfill -0.116843\\ \hfill -0.012096& \hfill 1.016352& \hfill 0.012924& \hfill 0.045527\\ \hfill 0.032553& \hfill 0.004589& \hfill 0.998258& \hfill 0.055794\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$,

$H_{tm2}^{-1}=\left[\begin{array}{rrrr}\hfill 1.012918& \hfill 0.166930& \hfill -0.011245& \hfill -0.209060\\ \hfill 0.071820& \hfill 1.028852& \hfill 0.106682& \hfill 0.048913\\ \hfill 0.065354& \hfill -0.028643& \hfill 0.978727& \hfill 0.171113\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$,

$H_{tm3}^{-1}=\left[\begin{array}{rrrr}\hfill 1.064280& \hfill 0.408684& \hfill -0.434687& \hfill 0.419931\\ \hfill 0.324798& \hfill 1.345943& \hfill 0.210470& \hfill 0.898905\\ \hfill 0.221952& \hfill 0.184071& \hfill 1.084612& \hfill -0.470323\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$,

$H_{tm4}^{-1}=\left[\begin{array}{rrrr}\hfill 0.112204& \hfill -0.176563& \hfill 0.677722& \hfill -1.436055\\ \hfill -0.037890& \hfill 0.631210& \hfill 0.318046& \hfill 0.523522\\ \hfill -0.659836& \hfill -0.011208& \hfill 0.253057& \hfill 0.588023\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$, and

$H_{tm5}^{-1}=\left[\begin{array}{rrrr}\hfill 0.599898& \hfill -0.757315& \hfill -0.454490& \hfill 1.325230\\ \hfill 0.730332& \hfill 0.576859& \hfill -0.516980& \hfill -0.534417\\ \hfill 0.809368& \hfill -0.039595& \hfill 0.673142& \hfill 2.244765\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$. The best of intermediate particles is $H_{tm1}^{-1}$ and it is better than $H_g^{-1}$, hence $H_g^{-1}=H_{tm1}^{-1}$. Then ${\mathit{x}}_{select}^b$ from Equation (20) is $H_{b\_select}^{-1}=\left[\begin{array}{rrrr}\hfill 0.995691& \hfill 0.122637& \hfill -0.006306& \hfill -0.116511\\ \hfill -0.012102& \hfill 1.016839& \hfill 0.012930& \hfill 0.045549\\ \hfill 0.030081& \hfill 0.212238& \hfill 1.000898& \hfill 0.065096\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$. After comparing $H_g^{-1}$ with $H_{b\_select}^{-1}$, $H_{b\_select}^{-1}$ is better, hence $H_g^{-1}=H_{b\_select}^{-1}$. Figure 3 depicts the registration mean-squared error (MSE) of the global best of the SR-PSO algorithm. We can see that, as the iteration proceeds, the MSE decreases and moves towards 3.22 × 10⁻³¹ at the 2000th iteration. Figure 4 shows the registration progress image of the global best particle in the 1st, 50th, 500th, and 2000th iterations. It can be seen that the registration result in the last iteration is an almost perfect match.

The registration mean-squared errors (MSEs) using the SR-PSO algorithm with and without refining with the ICP method are shown in

Table 3. The registration mean-squared error (MSE) in pixels² (the best value is in bold).

	α
	0.5	1.0	1.5	2.0
SR-PSO without ICP	7.68 × 10−2	9.23 × 10−2	3.22 × 10−31	9.68 × 10−2
SR-PSO with ICP	7.68 × 10−2	9.23 × 10−2	3.22 × 10−31	9.68 × 10−2

. The best result is when α = 1.5 with an extremely tiny error of 3.22 × 10⁻³¹. The best final registration result is also shown in Figure 5.

The best final transformation matrix is

$$H^{-1}=\left[\begin{array}{rrrr}\hfill 0.901895& \hfill 0.421703& \hfill 0.439339& \hfill 0.1\\ \hfill -0.143742& \hfill 0.826257& \hfill -0.237726& \hfill -0.1\\ \hfill -0.294063& \hfill 0.244980& \hfill 0.819157& \hfill 0.1\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right]$$

(24)

Although the MSEs of the final registration with the other values of α are not as good, the final registration images are almost aligned as shown in Figure 6. The reason for this mistake might be because the SR-PSO algorithm does not find the correct global minimum.

Now, we are ready to implement this algorithm on tooth reconstruction. We used the EinScan series, a commercial 3-D scanner, to collect a point-cloud data set in several consecutive views from two commercial tooth models. Examples of data acquisition setup and the point-cloud data are shown in Figure 7.

We implemented this system on two different tooth models, i.e., a regular tooth model and an orthodontic tooth model. We used only six consecutive point-cloud coordinate (x, y, z) views with an interval of 30 degrees in each model. The information of the tooth point-cloud data is shown in

Table 4. Tooth dataset information.

Model	Object View	Object Name	Points Number
Regular tooth model	1	Scan_0.asc	28,807
	2	Scan_1.asc	28,970
	3	Scan_2.asc	28,983
	4	Scan_3.asc	25,809
	5	Scan_4.asc	17,303
	6	Scan_5.asc	21,739
	Total	Six views	151,592
Orthodontic tooth model	1	Scan_0.asc	25,301
	2	Scan_1.asc	25,772
	3	Scan_2.asc	22,432
	4	Scan_3.asc	17,167
	5	Scan_4.asc	22,537
	6	Scan_5.asc	24,148
	Total	Six views	137,357

.

In this dataset, we randomly sampled each view into 60% of the original view. We also assumed that there was overlapping between each consecutive view. We also selected representative points inside the overlapping area using the voxel hull method [33,34,35] before we implemented the registration process with the parameter setting shown in Table 2. Since there were six consecutive views, we utilized the hierarchical registration to increase the registration performance as shown in Figure 8 with F = 6. At each level, we selected the best final registration result (SR-PSO algorithm with the ICP method) to survive to the next level of the hierarchical registration.

We also compared the results with those from the regular particle swarm optimization (PSO) algorithm with the ICP method. The registration process from the PSO algorithm with the ICP method is the same as the SR-PSO algorithm with the ICP method. However, the inertia weight update equation in PSO [29] is

$$w(t+1)=w_{\max }-t\frac{\left(w_{\max }-w_{\min }\right)}{\left(T\right)}$$

(25)

The parameter setting was similar to those used in the SR-PSO algorithm except that w_min, w_max, c₁, and c₂ were set to 0.4, 0.9, 2, and 2, respectively.

For the regular-tooth model, the registration MSEs for the SR-PSO algorithm and the PSO algorithm for each registration pair are shown in

Table 5. MSE from SRPSO and PSO for the regular-tooth model (the best value is in bold and the worst value is underlined).

	MSE in Micrometer2
	SR-PSO				PSO
View pairs	α = 0.5	α = 1.0	α = 1.5	α = 2.0
1 vs. 2	5.9300	6.1910	6.3274	6.3470	7.0328
2 vs. 3	5.0026	4.9070	4.8937	4.9507	5.3218
3 vs. 4	5.8824	5.4310	6.0370	5.5086	8.1904
4 vs. 5	5.2666	5.3131	5.2807	32.312	6.7879
5 vs. 6	6.0431	5.8163	5.9118	5.8166	6.5468

. Meanwhile, the registration MSEs for each of the two algorithms with the ICP method for each registration pair are shown in

Table 6. MSE from SR−PSO with iterative closet point (ICP) and PSO with ICP for the regular-tooth model (the best value is in bold).

	MSE in Micrometer2
	SR-PSO with ICP				PSO with ICP
View pairs	α = 0.5	α = 1.0	α = 1.5	α = 2.0
1 vs. 2	5.8628	5.8636	5.8632	5.8809	5.8156
2 vs. 3	4.8906	4.8912	4.8860	4.8884	4.8844
3 vs. 4	5.4030	5.4024	5.4017	5.4031	5.4055
4 vs. 5	5.1326	5.1253	5.1747	5.1261	5.1261
5 vs. 6	5.6829	5.6880	5.6828	5.6882	5.6956

. Figure 9 and Figure 10 show the best registration image of each consecutive pair from the SR-PSO algorithm without and with the ICP method, respectively. Again, the local optima are the cause of high MSE in pair-wise registration. Although Table 6 shows that the PSO algorithm with the ICP method is a little bit better than the SR-PSO algorithm with the ICP method in some pair-wise registrations, i.e., 1 vs. 2, 2 vs. 3, and 3 vs. 4, the difference is only in the order of 0.002 to 0.05 micrometer². When considering the final registration results of six consecutive views shown in

Table 7. MSE of the final registration of six consecutive views (micrometer²) for the regular-tooth model (the best value is in bold).

SR-PSO with ICP				PSO with ICP
α = 0.5	α = 1.0	α = 1.5	α = 2.0
7.3670	7.3668	7.3666	7.3667	17.1150

, we find that the best result from the SR-PSO algorithm with the ICP method and α = 1.5 is 7.3666 micrometer² whereas that from the PSO algorithm with the ICP method is 17.1150 micrometer². The result from the SR-PSO algorithm with the ICP method is better than that from the PSO algorithm with the ICP method. The final registration result of the regular-tooth model is shown in Figure 11. We can see that the final registration image provides a good visualization to human eyes.

For the orthodontic-tooth model, the registration MSEs for the SR-PSO algorithm and the PSO algorithm for each registration pair are shown in

Table 8. MSE from SR-PSO and PSO for the orthodontic-tooth model (the best value is in bold).

	MSE in Micrometer2
	SR-PSO				PSO
View pairs	α = 0.5	α = 1.0	α = 1.5	α = 2.0
1 vs. 2	5.5553	5.7212	5.6724	5.7199	28.973
2 vs. 3	6.2493	6.2868	6.1613	6.2224	35.490
3 vs. 4	179.95	12.719	5.4687	12.779	6.7116
4 vs. 5	7.0814	6.7577	6.5847	8.6855	14.036
5 vs. 6	5.4288	5.3262	5.3638	5.4472	35.000

. Figure 12 shows the best registration images of different pairs from the SR-PSO algorithm. Meanwhile, the registration MSEs for the SR-PSO algorithm with the ICP method and the PSO algorithm with the ICP method for each registration pair are shown in

Table 9. MSE from SR-PSO with ICP and PSO with ICP for the orthodontic-tooth model (the best value is in bold).

	MSE in Micrometer2
	SR-PSO with ICP				PSO with ICP
View pairs	α = 0.5	α = 1.0	α = 1.5	α = 2.0
1 vs. 2	5.5093	5.5101	5.5093	5.5127	5.5130
2 vs. 3	6.1441	6.1440	6.1440	6.1444	6.1570
3 vs. 4	5.2854	12.511	5.2706	12.514	5.2847
4 vs. 5	6.3945	6.3948	6.3946	6.3999	6.3933
5 vs. 6	5.2859	5.2801	5.2810	5.2814	5.2815

. The best registration images of different pairs from the SR-PSO algorithm with the ICP method are shown in Figure 13. We can see that using only the SR-PSO algorithm provides better results than using only the PSO algorithm. However, the PSO algorithm with the ICP method is better than the SR-PSO algorithm with the ICP method in one of the pair-wise registrations (4 vs. 5) with a little difference. The final registration MSEs of the six consecutive views are shown in

Table 10. MSE of the final registration of six consecutive views (micrometer²) for the orthodontic-tooth model (the best value is in bold).

SR-PSO with ICP				PSO with ICP
α = 0.5	α = 1.0	α = 1.5	α = 2.0
7.4130	7.4141	7.4141	7.4131	7.4672

. The best final 3-D registration image is from the SR-PSO algorithm with the ICP method and α = 0.5. It has an MSE of 7.4130 micrometer² whereas the comparable error of 7.4672 micrometer² is achieved by the PSO algorithm with the ICP method. We can also see in Figure 14 that the final 3-D registration image of the orthodontic-tooth model also provides a good visualization to human eyes.

4. Conclusions

To aid in orthodontics or the diagnostic and treatment planning process in dental caries, a 3-D tooth model reconstruction is one of the important parts. In this paper, we introduce a registration process using our proposed statistical randomization-based particle swarm optimization (SR-PSO) algorithm with the iterative closet point (ICP) method to determine the optimal affine transform between images in 3-D. Because there were several consecutive images scanned from a commercial tooth model, we implemented a hierarchical registration process to produce the final registration image. We also compared the results with those from the particle swarm optimization (PSO) method with the ICP method. The best results from the SR-PSO algorithm with the ICP method for the regular-tooth model and the orthodontic-tooth model yielded the mean-squared error (MSE) of 7.3666 micrometer² and 7.4130 micrometer², respectively. For the sake of comparison, the MSEs from the PSO algorithm with the ICP method for the two tooth models were 17.115 micrometer² and 7.4672 micrometer², respectively. The SR-PSO algorithm with the ICP method provided better results than that of the PSO algorithm with the ICP method. Moreover, the resulting 3-D images from the SR-PSO algorithm with the ICP method are also viewable by human eyes and are useful for experts. In future work, we will implement this algorithm for the tooth models with defections in order to simulate dental caries in real situations.