Abstract
Multimodality image registration and fusion has complementary significance for guiding dental implant surgery. As the needs of the different resolution image registration, we develop an improved Iterative Closest Point (ICP) algorithm that focuses on the registration of Cone Beam Computed Tomography (CT) image and high-resolution Blue-light scanner image. The proposed algorithm includes two major phases, coarse and precise registration. Firstly, for reducing the matching interference of human subjective factors, we extract feature points based on curvature characteristics and use the improved three point’s translational transformation method to realize coarse registration. Then, the feature point set and reference point set, obtained by the initial registered transformation, are processed in the precise registration step. Even with the unsatisfactory initial values, this two steps registration method can guarantee the global convergence and the convergence precision. Experimental results demonstrate that the method has successfully realized the registration of the Cone Beam CT dental model and the blue-ray scanner model with higher accuracy. So the method could provide researching foundation for the relevant software development in terms of the registration of multi-modality medical data.
1. Introduction
With the rapid development of the image processing, reverse engineering, computer-aided design, and image guidance play an important role in dental implant surgery [1]. A number of dental and jaw imaging modalities, including video imaging [2], Computed Tomography (CT) [3], and magnetic radiotherapy imaging (MRI) [4], are used in image-aided dental implant surgery. Performing the multimodality image registration from different imaging devices has complementary significance for guiding dental implant surgery. When dentists design schemes, they want to consider the morphology of the soft tissue and the hard tissue together. CT image and visible spectrum image fusion models [5] can multiply various factors and help the dentists achieve the optimal implant scheme.
Traditional registrations of Cone Beam Computed Tomography (CBCT) and other three-dimensional scanner models mainly rely on the fiducial markers method, which need to manually select the feature points. In [6–8], manual feature selecting methods have been tried for the registration of the mandible CT and plaster cast of the dental models. Because the images are complex and the location accuracy is limited by the experience of the operator and the operating state, the manual feature selecting method might provide misleading information and lead to instability problems. Noh et al. [9] designed an occlusal form board, with titanium spheres markers, to match the image using the titanium spheres markers as feature points. In this model, the position of makers is stable, but the drawback is the extra burden for patients. The methods mentioned above cannot perfectly meet the needs of developing image-aided dental surgery. Some attempts have been made to select the feature point set automatically. Two main difficult problems are taken out in the process of the registration; one is the feature point selection and the other is the merge of different modality images with different resolutions.
Iterative Closest Point (ICP), in recent years, many researchers have extended or improved the traditional ICP algorithm. Xie et al. [10] proposed the fractional ICP (FICP) algorithm to trim abnormal points effectively, which enhanced the robustness of the ICP algorithm. Choi et al. [11] proposed the improved - trees traversal method to accelerate finding process of the nearest points. Bae and Lichti [12] further extended the improved ICP method based on the boundary feature points of the point cloud, which improved the efficiency and accuracy of traditional ICP algorithm greatly.
The aim of this work is to develop an improved Iterative Closest Point (ICP) algorithm in terms of the problem of registration with different image resolutions. Firstly, extract feature point set based on curvature characteristics, and then achieve the registration results by two procedures: the coarse registration and the precise registration. To verify the reliability of the improved ICP algorithm, we test it with 3D CBCT dental model and inEos Blue-light scanned model. These two models provide reconstructed CT dental data and visible light scanned dental-surface data, respectively. Experimental results show that the accuracy is increased significantly and this solution meets the requirements of the image-assisted surgery.
The rest of this paper is organized as follows. Section 2 provides an overview of the ICP algorithm and proposes an improved ICP algorithm. Section 3 describes the application of the improved ICP algorithm. Section 4 reports the results of experiments conducted on the CBCT model and the inEos blue scanned model and provides the discussions. Conclusions are presented in Section 5.
2. The Improved ICP Algorithm
Generally, different medical imaging equipment differs in imaging mechanism, which contributes to different emphases on image information and makes image resolution difference. Thus, it results in complex difficulties in the image fusion; for example, the preprocessed 3D CBCT dental model contains two parts: dental crowns and roots, while the 3D inEos Blue-light dental model contains dental crowns only. Additionally, the resolution of the CBCT model is 0.625 mm, while that of inEos Blue-light scanner is 0.1 μm. In view of data characteristics above, an improved ICP algorithm is proposed in this paper. First of all, feature point set of the three-dimensional surface images was extracted based on curvature characteristics to realize the coarse registration and to achieve the transformation matrix. Then, the feature point set and reference point set, obtained by the initial registered transformation, are processed in the precise registration step using ICP.
2.1. ICP Algorithm
ICP algorithm is commonly used as registration method, which enables the minimum mean square error in the distance of corresponding point or point surface between two point sets approaching the minimal by iterative calculations [13]. It repeats the process of calculating the optimal rigid transformation of the corresponding point set until convergence criterion is satisfied; that is, it has met the precision of registration. The transformation relation is shown as function (1), and the convergence criterion is shown as function (2): where and present point sets corresponding to two images, respectively, is a rotation matrix, is a translation vector, and is the error. When the error reaches the minimum, the convergence criterion meets requirement [14].
The traditional ICP algorithm used in registration requires one of the two point sets to be the subset of the other one. The relatively initial positions of two point sets are limited strictly in ICP, which ensure no big gap between the two point sets. When the constraints are unsatisfied or the huge deviation exists, it will affect the processing results and make the registration unreliable.
2.2. The Selection of Feature Point Sets
The data of the CBCT dental model and the inEos Blue scanned model are both stored in STL file form and 3D models. The image surface discrete curvatures are sampled as the feature set. Because the STL Model is comprised of many triangular meshes, an improved Voronoi method is implemented to calculate vertex curvature [15].
As shown in Figure 1, is a vertex of the network model. The curvature of can be calculated by where is the angle between the side and the side , and is the area of all triangle meshes which are connected with .
Set as the threshold. When the curvature of one point is greater than , it will be marked as feature point and put it into alternative feature list. The feature point set is defined as .
2.3. The Coarse Registration
The improved three point’s translational transformation method is used in this paper based on curvature characteristics to calculate rotation matrix and translation vector , so as to realize the semiautomatic extraction of corresponding points for two dental image models. At the same time, this method can reduce the interference of human factors in the matching process and enhance the accuracy of coarse registration.
2.3.1. The Improved Three Point’s Translational Transformation Method
For the two image models, one is defined as reference model and the other is defined as target model. Select three points, whose curvatures in the top 10 percent and distributions are homogeneous, as the reference points from the reference model and the corresponding points are in target model. After achieving reference points, we calculate the curvature of each vertex in target model and compare it with their value of the reference points. When the error is smaller than the setting threshold, mark this point as the centralized point of corresponding points and put it into corresponding points list. Due to the fact that there are some feature points with the same or similar curvature, it is common to lead to one reference point corresponding to some points. So we construct more powerful constraint condition, the distance of three points, and the area of triangle of three points, further to select the corresponding point set. This constraint condition is independent of rigid transformation and will improve the reliability and efficiency of determining corresponding points, and its specific steps are as follows.
Assume there are three known reference points , , and and three corresponding points , , and that need to be determined. Let , , and represent the distance between , or , and , and , respectively. is the area of triangle composed of , , . The computing function of the distance and the area are shown as follows: The distance constraint condition is: where and are the distances between any two points of the reference points , , and , and the corresponding points , , , respectively. is a threshold defined. The area constraint condition is shown by: where is the area of triangle composed of , , and ; is the area of triangle composed of the three corresponding points , , and ; and is the set error.
Through the above stated constraint conditions, we can select the corresponding points , , and of the reference points , , and .
2.3.2. The Course Registration
Local coordinate system of one of the two image models is built by reference points in the reference model. Assume that is the origin of coordinate and the direction from to is axis direction. Local coordinates can be expressed by , where
Similarly, the local coordinate built by corresponding points , , and is . By formula deducing, we could obtain rotation matrix and translation vector corresponding to the transformation of coordinate to :
Then, the matrices and are applied to the two target point sets and we could realize the course registration.
2.4. The Precise Registration
The accuracy of coarse registration cannot satisfy the ultimate requirements [16, 17], so ICP algorithm is adopted for getting higher registration accuracy.
Assume that the feature point sets of Blue-light scanner dental model and CBCT dental model are and , respectively, where and stand for the number of the elements . , are the rotation matrix and translation vector, which are obtained through the coarse registration. Assume that and stand for the iterate times and the precision threshold. Treat as the target point set and as the reference point set. Then carry out the ICP algorithm to make precise registration that consists of the following stages.
Step 1. Given that , , , , and , use - tree method to calculate the corresponding point set , which belongs to , and each point is with the least Euclidean distance to the point . All the corresponding points must be unique. Define the initial residual error; that is, .
Step 2. Use and transform point set to obtain the point set , where . Then calculate the corresponding point set , which also belongs to , and each point is with the least Euclidean distance to the point . Then compute the residual error .
Step 3. Judge whether the precision meets the requirements after iterating times. When , is satisfied stop iterating. Otherwise, continue to Step 4.
Step 4. Calculate the rotation matrix and translation vector of transform relation between and through the Singular Value Decomposition (SVD) method [18]. Then renew , and transfer to Step 2.
The Singular Value Decomposition (SVD) method, used to calculate the rotation matrix and translation vector, is demonstrated in the appendix.
Figure 2 is the flow chart of precise registration.
3. The Registration of CBCT Dental Model and inEos Blue Dental Model
The experimental procedures include acquiring data, preprocessing, extracting feature points, and registration.
3.1. Data Acquisition and Preprocessing
The main procedures of data acquisition and preprocessing are shown in Figure 3. Firstly, using 3D Sirona inEos Blue Blu-ray scanner (resolution is 0.02 mm) patient’s plaster dental model is scanned. Next, the own software of scanner is applied to obtain STL format dental file. Geomagic software is used to preprocess the exported file, and then, export the preprocessed file, which is STL that format includes dental crown only. The result is shown in Figure 4.
Sirona GALILEOS system (resolution is 0.1 mm) is adopted to scan the same patient’s mandible to obtain the CT scanning data. The data is DICOM format and then it is imported into commercial processing software called Mimics, which is used to reconstruct 3D CT model. After preprocessing, the 3D mandible model with STL format is exported. Figure 5 shows the CBCT dental model with both dental crown and root.
3.2. Extracting the Feature Points
The inEos Blue dental model, showed in Figure 4 as the sample, set the threshold . Calculate the vertex curvature of the dental crown using the improved Voronoi method mentioned before. When the curvature value of any point is more than , mark it as a feature point and put it into point set list, which is denoted by . The red dots showed in Figure 6 represent the feature points.
3.3. Coarse and Precise Registration
The resolution of 3D inEos Blue dental data is higher than CBCT data, while CBCT model includes the integral dental data. The improved three point’s translational transformation method is used to realize coarse registration, and the results are shown in Figure 7. Select three corresponding points , , and as reference points from CBCT dental model and set the threshold 0.00015. Then extract the corresponding points , , and and calculate rotation matrix and translation vector . The results are listed in Table 1.
Formula (8) is used to realize coarse registration, and Figure 8 is the result of coarse registration.
Based on the results of coarse registration, the precise registration is carried out for improving the accuracy of registration further. Figure 9 is the result, which reveals that the two models can be fused more accurately.
4. Discussions
Supported by the VC6.0 platform, we used VC++ and OpenGL to conduct simulation and experiments.
4.1. The Accuracy of Registration
We show dental models registration results to evaluate the accuracy of our algorithm, and the average distance of corresponding points is calculated. For the first experiment the average distance of the coarse registration, shown in Table 2, is 0.7291 mm while it is reduced to 0.3835 mm after precise registration.
Other four groups of patient’s dental image data are carried out by the same procedure, which also conclude inEos Blue scanner and CBCT 3D dental models and yield the same conclusion. The total average distance of the precise registration is 0.8266 mm, compared with 1.4486 mm of the coarse registration. Comparing with the coarse registration accuracy, the precise one is improved significantly, which basically meets the requirements of helping dentist for designing remedies.
4.2. The Reliability and Advantages of the Feature Extraction
Considering the characteristic of the STL model data itself, we take the CBCT 3D dental model as the reference model and adopt the constraint-based curvature method. As mentioned earlier, extract the corresponding feature points automatically in the inEos Blue scanner dental model. Experimental results show that the method could extract the corresponding feature points correctly and be operated simply and reliably. Moreover, it efficiently overcomes the shortages of both the traditional surface-based marking method and the corresponding manually selected feature points method.
4.3. Further Extension and Improvement
Multimodality image registration and fusion technology can make full use of the complementary characteristics which are provided by different imaging models. The technology could perfectly integrate the image data to provide abundant information for clinical diagnosis. The improved ICP method proposed in this paper realized the feature points extraction and image registration successfully between the inEos Blue scanner dental model and the CBCT 3D dental model. Besides, it provides reliable registration method in the computer-aided dental implant design system. The research and the method in this paper will be the reference and application for similar studies and other multimodality image registrations.
Imperfectly, in the process of extracting feature points with curvature characteristics, the threshold has to be appropriately adjusted according to the model complexity and the accuracy requirements, which make the registration accuracy to a larger degree depend on the reasonable choice. The adaptive threshold selection method demands our further research.
5. Conclusion
As the traditional ICP algorithm cannot realize the image registration with different resolution, we propose an improved ICP algorithm in this paper. The procedure of this algorithm used in the registration of 3D CBCT dental model and 3D inEos Blue scanning dental model is discussed in detail. The experimental results verified that this method can guarantee a high accuracy in registration.
Appendix
Assume the unit quaternion is . The quaternion expression after coordinate system unit vector rotated angle is The corresponding rotation matrix is
Assume the centroids of two point sets and are and , respectively. The covariance matrix between the point set and is , and represents the trace of matrix . Define symmetric matrix :
Using Jacobi method, calculate the eigenvalue and the eigenvector of matrix . The corresponding eigenvector of the maximum eigenvalue is the quaternion , which meets the requirement. As these are known, then we can calculate rotation matrix , that is, ; and the translation vector is .
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank Professor Jianping Geng for providing some of the dental images. This work was supported by a Grant from the Qing-Lan Project of Jiangsu Province.