Computerized Medical Imaging and Graphics 44 (2015) 62–68
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Automated 3D–2D registration of X-ray microcomputed tomography with histological sections for dental implants in bone using chamfer matching and simulated annealing Kathrin Becker a,b,∗ , Martin Stauber c , Frank Schwarz b , Tim Beißbarth a a
Department of Medical Statistics, Biostatistics Group, University Medical Center, Georg-August University, Humboldt Allee 32, 37073 Göttingen, Germany Department of Oral Surgery, Westdeutsche Kieferklinik, Heinrich-Heine University, Moorenstr. 5, 40225 Düsseldorf, Germany c Scanco Medical AG, Fabrikweg 2, 8306 Brüttisellen, Switzerland b
a r t i c l e
i n f o
Article history: Received 15 November 2014 Received in revised form 1 March 2015 Accepted 17 April 2015 JEL classification: I10 C65 Keywords: Image registration Micro-computed tomography Histology Osseointegration Dental implants Histomorphometry
a b s t r a c t We propose a novel 3D–2D registration approach for micro-computed tomography (CT) and histology (HI), constructed for dental implant biopsies, that finds the position and normal vector of the oblique slice from CT that corresponds to HI. During image pre-processing, the implants and the bone tissue are segmented using a combination of thresholding, morphological filters and component labeling. After this, chamfer matching is employed to register the implant edges and fine registration of the bone tissues is achieved using simulated annealing. The method was tested on n = 10 biopsies, obtained at 20 weeks after non-submerged healing in the canine mandible. The specimens were scanned with CT 100 and processed for hard tissue sectioning. After registration, we assessed the agreement of bone to implant contact (BIC) using automated and manual measurements. Statistical analysis was conducted to test the agreement of the BIC measurements in the registered samples. Registration was successful for all specimens and agreement of the respective binary images was high (median: 0.90, 1.–3. Qu.: 0.89–0.91). Direct comparison of BIC yielded that automated (median 0.82, 1.–3. Qu.: 0.75–0.85) and manual (median 0.61, 1.–3. Qu.: 0.52–0.67) measures from CT were significant positively correlated with HI (median 0.65, 1.–3. Qu.: 0.59–0.72) between CT and HI groups (manual: R2 = 0.87, automated: R2 = 0.75, p < 0.001). The results show that this method yields promising results and that CT may become a valid alternative to assess osseointegration in three dimensions. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Biological processes associated with the osseointegration of endosseous dental implants have been extensively investigated in numerous preclinical and clinical studies [1–3]. The level of osseointegration is typically evaluated by histological assessment of bone to implant contact (BIC) areas in undecalcified hard tissue sections [4]. One major drawback of nondecalcified sectioning, however, is a significant loss of biological
∗ Corresponding author at: Department of Orthodontics, Westdeutsche Kieferklinik, Heinrich-Heine University, Moorenstr. 5, 40225 Düsseldorf, Germany. Tel.: +49 211 8118140; fax: +49 211 8116550. E-mail addresses: [email protected] (K. Becker), [email protected] (M. Stauber), [email protected] (F. Schwarz), [email protected] (T. Beißbarth). http://dx.doi.org/10.1016/j.compmedimag.2015.04.005 0895-6111/© 2015 Elsevier Ltd. All rights reserved.
tissues during the histological (HI) processing of serial sections, and therefore, only a maximum of 2–3 sections are available which may not reflect the entire healing processes at the investigated site. Another drawback is the amount of time needed for the processing of hard tissue sections, delaying data acquisition for several months after completion of the experimental procedures. Moreover, staining of specimens is a highly sensitive process and commonly associated with marked variances in color intensity. Furthermore, evaluation is limited to 2D, so different orientations can not be investigated, and killing of the animal is needed to harvest biopsies for histology thus for each time point different animals are needed. In recent years, micro computed tomography (CT) has become a promising tool for quantitative three-dimensional (3D) analyses of new bone formation following implant placement and various types of bone augmentation procedures. This non-destructive and non-invasive volumetric analysis offers many advantages including
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
a fast and automated characterization of structural features such as bone volume density [5,6] or the assessment of BIC values [7,8] in different planes and even in 3D. Furthermore, with the advent of intravital CT scanners, different time points can be studied within the same animal. Despite these advantages, assessment of BIC values in CT scans may be impeded by beam hardening artifacts, partial volume effects and scattering [9,10] which occur in the proximity of dental titanium implants. While partial volume effects can be eliminated by performing BIC measurements few voxels away from the implant, beam hardening and scatter artifacts may distort BIC measurements. Due to the lacking proof of comparability, CT analyses are usually complemented with conventional HI at the moment [11–13]. To verify if CT and HI yield the same results, and to validate novel BIC measurements in the future, 2D corresponding slices have to be compared with the gold standard HI. For this a 3D–2D registration is required that aligns corresponding image elements from CT and HI. Few approaches for 3D–2D registration are described in literature, mostly reporting on the registration of CT or MRI volumes with projected 2D X-ray [14,15], ultrasound [16,17] or fluoroscopic images [18–20]. Registration of single 2D slices obtained during interventional procedures with CT [21] or HI section with CT did not gain much attention so far since the registration relies on manyto-one correspondence of the underlying 3D-to-2D mapping as well as a statistical consistency problem owing to the reliance of a similarity measure on small sampling size of an input 2D image [22] and big differences in information content [23]. Manual approaches for 3D–2D registration are hard to realize because HI is usually not located parallel to a reference plane within CT, thus, a seven dimensional optimization problem needs to be solved (rotation, translation and scaling). Sarve et al. were the first to report on automated computerized registration for synchroton radiation microcomputed tomography and histology using the model of osseointegrated orthopedic titanium (TI) implants [10]. They used normalized mutual information to assess similarity and assumed HI being located directly at the implant axis. They also tried to employ simulated annealing without a prior first guess and fast cooling so a failure rate of 25% was reported by the authors [10]. The aim of the present study was to develop a robust automated computerized registration method for HI images with CT volumes. The specific goals were (i) to establish an algorithm to automatically register 2D HI images with a 3D CT volume, and (ii) to assess the agreement of BIC values in the registered images.
2. Acquisition of data 2.1. Tissue samples The specimens (n = 10) used for the present analysis were randomly selected from tissue biopsies obtained in a recently published experimental study [24]. In brief, each biopsy contained one commercially available titanium implant (Ø 3.8 mm, length: 9 mm, Conelog® Screw-Line Implant, Camlog Biotechnologies AG, Basel, Switzerland) connected with a diameter-reduced wide body Ti (Ti6Al4V) healing abutment (Ø 3.8 mm, height: 4 mm, Camlog). These biopsies were retrieved at 20 weeks after implant placement and non-submerged healing in the mandible of foxhounds. The study protocol was approved by the appropriate local authority (Landesamt fr Natur und Verbraucherschutz, Recklinghausen, Germany). The biopsies were dehydrated using ascending grades of alcohol and xylene, infiltrated and embedded in methylmethacrylate (Technovit 9100 NEU; Heraeus Kulzer, Wehrheim, Germany),
63
which is the standardized procedure for non-decalcified sectioning. This procedure takes 20 h, during which negative influence of polymerization heat was avoided by storing the specimen at cold atmosphere (4 ◦ C). 2.2. MicroCT acquisition After embedding in methylmethacrylate, the biopsies were scanned with a CT 100 (Scanco Medical AG, Brttisellen, Switzerland). The X-ray tube was operated at an energy of 90 kVp and an intensity of 88 A. The integration time was chosen to 1200 ms and the frame averaging was set at 2×. The data was reconstructed into 3D volumes with an isotropic nominal resolution of 8.6 m. 2.3. Histological preparation and image acquisition The technique has been described in detail [24]. Briefly, after CT analysis, histological preparation was performed. Each implant site was cut in the vestobulo-oral direction along with the long axis of the implant using a diamond band saw (Exact(r); Apparatebau, Nordersted, Germany). Serial sections of approximately 300 m in thickness were prepared, glued with acrylic cement (Technovit 7210 VLC; Hareus Kulzer), to silanized glass slides (Super Frost; Menzel GmbH, Braunschweig, Germany) and ground to a final thickness of approximately 40 m. Subsequently, from each tissue block, one specimen was randomly selected for the present analysis and stained with paragon. For image acquisition a color CCD camera (Color View III, Olympus, Hamburg, Germany) mounted on a binocular light microscope (Olympus BX50, Olympus, Hamburg, Germany) was used. Digital images for registration (original magnification 100×) were obtained using the Cell D program (Cell D® , Soft Imaging System, Münster, Germany). Additionally, the specimens were scanned with Aperio ScanScope Slide Scanner (Aperio Technologies Inc., CA, USA) for BIC assessment (original magnification 200×) and the region of interest was exported as snapshot using Aperio ImageScope v11.1.2.752 (Aperio Technologies). 3. Image registration 3.1. Image preprocessing and feature extraction The image preprocessing and registration was performed with Matlab (Matlab R2013a 64-bit, The Mathworks Inc., Massachusetts, USA). During preprocessing, both HI and CT were each transformed into two binary images showing either the dental implant or the bone (Fig. 1). To segment the implant from HI, thresholding and component labeling was applied. The edges were detected using Canny edge filter. For segmentation of the bone tissue from HI, a combination of thresholding, component labeling and morphological filters (erosion, dilation) with a 3 × 3 disk as structuring element was used. On the CT slices, the thresholds were chosen manually such that either the bone tissue or dental implant was segmented. 3.2. Adaption of image resolution (scaling) The element size of histology was obtained from a bar of fixed length (5 mm) located within the histological images. The CT element size was read from the header of the CT-files given in Scanco Medical AIM-format. Adaption was achieved by scaling HI images with the ratio of both element sizes, i.e. sizehi /sizect .
64
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Fig. 1. Image preprocessing procedure: to segment the bone tissue from HI, a mask to cover the bone was computed using thresholding and morphological filters. With component labeling, the marrow spaces could be cleared. Segmentation of the implant from HI was conducted by finding the largest component. At the CT images segmentation was done using thresholding.
the approximated Euclidean distance to the nearest edge, the alignment score was given by the root mean squared average (RMS), which averages the intensities xi to which the ei were mapped:
n 11 RMS = xi n
3
(1)
i=1
The division by 3 is required since the 3–4 Chamfer distance was used instead of the Euclidean distance. The optimal match was found when no transformations were possible that would reduce RMS. Convergence of CM was shown by Borgefors [25]. Due to proximate rotation symmetry of the implant, RMS from CM with slices located around the implant axis was expected to be likewise. Thus the segmented bone from HI was mapped to the segmented bone from CT slice using the transformation matrix from CM. Then, similarity L was measured between the binary images by averaging the L1 norm and subtraction from 1 (xi , yi denote the values of aligned pixels): Fig. 2. Chamfer matching principle: the distance transform (DT) of 72 slices, that had been extracted from the 3D volume before and were located around the implant axis, was computed. The output was the optimal alignment of the corresponding 2D edges.
3.3. First guess registration using chamfer matching and principal component analysis Chamfer matching (CM) was used to find a first guess for the start of the fine registration and was implemented as described by Borgefors [25] for the 2D case (Fig. 2). Two pyramid levels were used and initialization was achieved by finding the implant axis from CT, which was approximated with the first axis (PA) from principal component analysis (PCA): to conduct PCA, all edge vectors from HI and CT were centered and stored as row vectors vi in a matrix. Then, its covariance matrix and eigenvalue decomposition were computed. The eigenvector with largest eigenvalue yielded PA. Next, a vector orthogonal to the implant axis (third principal axis, i.e. eigenvector with smallest eigenvalue) was rotated in 5 degree steps around PA and a series of 72 slices was extracted from the CT volume. From each slice, the 3–4 distance transform (DT) was computed [25] and the implant edges ei from HI were moved along DT until an optimal match was found. Since each intensity of DT reflects
1 xi − yi 1 n n
L =1−
(2)
i=1
Hence, RMS was used to find the transformation of HI and a slice extracted from CT. L was used to compare similarity of different slices extracted from CT with HI. These slices were obtained by rotating a vector orthogonal to the implant axis in 5 degree steps. At the position with highest L score, additional 1 degree rotations of the normal were performed (10◦ to each direction). Missing intensity values of oblique slices were interpolated using trilinear interpolation. The parameters revealing the optimal L score finally yielded the first guess registration of the implants using CM. 3.4. Fine registration of the segmented bone using simulated annealing Simulated annealing (SA) is a robust optimization method incorporated from thermodynamics. It was implemented with 6 degrees of freedom (rotation, translation) to find a slice from CT that maximizes L when being aligned with HI. Input were the segmented bone binary images from CT and HI as well as the so far optimal position and normal vector from Chamfer Matching. 3D–2D
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Fig. 3. Manual BIC measurements: a segmented line tool was used to measure the lengths of the implant outline (left) and the implant fractions being in direct contact with bone (right, green arrow) below the implant shoulders. The sum of the fractions divided by the length of the implant outline was defined as BIC measure. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
65
Fig. 4. Automated BIC measurements: after thresholding (right image), a mask was located 3 pixel distant from the implant outline. The number of pixels classified as bone (red pixels) divided by the entire number of pixels within the mask (red or green pixels) yielded the automated 2D BIC measure. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
matching was achieved by refining the position and normal vector of the CT slice, while HI remained fixed. The algorithm uses a parameter temperature T that was initialized with 300 and cooled down to 30. The step width for transformations was set to s = T/100. The Boltzmann constant k was set to 0.01. The energy states E1 , E2 were set to the alignment scores L of the present and next slice to be extracted from the volume. The next transformation was chosen randomly and excepted if the probability p of a system to change energy from E1 to E2 was smaller than a random number. p=
⎧
⎨ exp − E1 − E2 if E2 ≤ E1 ⎩
kT
1
if E2 > E1
The temperature was cooled down by 1% after 250 iterations and parameters yielding greatest L finally provided the optimal alignment. The code of the implementation can be found at: http:// kathrin-becker.eu/code/Registration.zip
Fig. 5. Threshold intervals were selected such that the implant as well as bone tissue could be segmented clearly from CT. However, in proximity to the implant (red arrow, right image) artifacts were sometimes mis-classified as bone (green arrow, left image) and thus increased automated BIC measures (colored mask, left image). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
3.5. Quantification of registration results and BIC measurements For quantification of the registration results, the L scores from registration were used. In addition, the registration results were visually examined and to evaluate agreement of BIC between the two imaging modalities, BIC was assessed in the region below the implant shoulder both automated and manually. This measurement refers to the relative fraction of bone being in direct contact with the implant. Manual measurements were performed using ImageJ (National Institutes of Health, USA) by two experienced investigators (histology: KB, CT: MS). The length of the implant contour and the lengths of the implant outline in direct contact with the bone were measured in the region below the implant shoulders using the Segmented Line measurement tool (Fig. 3). Prior to the start of the analyses, a calibration procedure was initiated for the investigators and revealed that repeated measurements were similar at >95% level. Automated measurements to obtain CT BIC were implemented and performed using Matlab R2013a. The automated BIC was assessed within a mask of 2 pixel width located below the implant shoulder and 3 pixel distant from the implant surface (Fig. 4). The mask was computed by applying image dilation and erosion filters
on the segmented implant. The position of the implant shoulders was specified manually using Matlab’s impixel feature. The BIC measure was given by the relative fraction of pixels classified as bone. To segment the bone, the same procedure as described in Section 3.1 was used. 3.6. Statistical analysis The statistical analyses were performed using the R statistics software version 2.12.2 [26]. Test for linear association between 2D HI BIC and 2D CT BIC was performed by calculation of linear models. Results were considered statistically significant at p < 0.05. 4. Results 4.1. Image preprocessing Segmentation of the implants from CT was successful for all specimens investigated. Segmentation of bone tissue from CT was sometimes impeded by artifacts, which were classified as bone in close proximity to the implant surface (Fig. 5).
66
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Table 1 The registration results are presented column-wise. Similarity scores L: chamfer matching (CM), simulated annealing (SA). BIC measurements: histology (HI), manual CT measurements (m), automated CT BIC measurements (a). The last rows shows the respective mean values ± standard deviation, median and quartiles. Specimen
L (CM)
L (SA)
HI BIC
CT BIC(m)
CT BIC(a)
1 2 3 4 5 6 7 8 9 10
0.86 0.83 0.83 0.88 0.83 0.87 0.83 0.80 0.87 0.82
0.91 0.90 0.93 0.93 0.90 0.91 0.85 0.87 0.89 0.92
0.60 0.77 0.73 0.80 0.46 0.64 0.65 0.59 0.22 0.72
0.61 0.79 0.68 0.75 0.52 0.61 0.55 0.50 0.26 0.62
0.86 0.82 0.86 0.91 0.75 0.89 0.77 0.74 0.61 0.81
Mean ± SD Median 1.–3. Qu.
0.85 ± 0.03 0.83 0.83–0.87
0.9 ± 0.03 0.91 0.9–0.92
0.62 ± 0.17 0.65 0.59–0.73
0.59 ± 0.15 0.61 0.52–0.67
0.8 ± 0.09 0.82 0.75–0.85
Fig. 7. The registration of specimen 7 was reliable even though parts of the bone tissue were damaged on production and despite visible color variance of the staining agent. However, this was reflected in the lowest L score of 85 %. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
4.3. Assessment of bone to implant contact Significant correlation of BIC for both, manual (Pearson = 0.94, p < 0.001, CI: 0.80–1) and automated measurements were found (Pearson = 0.88, p < 0.001, CI: 0.65–1). Regression analysis yielded a significant association between HI and CT groups (manual mea2 surements: regression coefficient ˇ = 1.01, Rkorr = 0.87, p < 0.001; 2 = automated measurements: regression coefficient ˇ = 1.7, Rkorr 0.75, p < 0.001). BIC values were higher for the automated measurements (median: 0.82; 1.–3. Qu.: 0.75–0.85) compared to the manual measurements (median 0.61; 1.–3. Qu.: 0.52*-0.67) and HI (median 0.65; 1.–3. Qu.: 0.59–0.72). Artifacts were present mainly in the apical portion of the CT slices. Fig. 6. Registration results of specimen 6: this specimen was not sawed directly at its center thus chamfer matching could only provide a rough approximation. Interestingly, the marrow space at the arrow is hardly visible in CT.
The implants could be segmented clearly from HI without adjustment of threshold levels. The bone tissue segmentation from HI required minor adjustments of threshold levels due to low staining intensities in three specimens. As small components were removed during segmentation, very small bone elements were missing after segmentation. It was found that this segmentation procedure was very sensitive to the staining process and may need to be adjusted for different staining procedures.
4.2. Image registration Computation of the implant axes and subsequent application of chamfer matching provided reliable first guesses when using similarity of the segmented bone L for inter-slice comparison. For the 2D matching of the segmented bone, RMS was in the range of 0.47–4.10. When comparing the RMS values obtained by rotating the normal around the implant axes, no clear global optimum could be found for RMS scores among different slices due to shape similarities of the implants around their axes. The transformation matrix from CM was used to match also the bone tissue which, in contrast to RMS, yielded clear optima and thus a reliant first guess for SA. The similarity scores L were on median 0.83 with 1st to 3rd quartile reaching from 0.83 to 0.87 (Table 1). Simulated annealing provided successful fine registration for all specimen and increased the similarity scores L in all cases (median 0.9, 1.–3. Qu.: 0.89–0.91) (Table 1). Visual examination corroborated the high agreement of marrow spaces and bone crest (Fig. 6).
5. Discussion In this study, we presented a novel method for automated 3D–2D registration of X-ray microcomputed tomography with histological sections. This method was developed for the direct comparison of these two techniques, and to evaluate the agreement of the bone to implant contact (BIC). To prove that the correct position and normal vector for the corresponding CT slice were found, the results were quantified by means of similarity scores and by visual examination. To the best of the authors’ knowledge, the only study dealing with registration of HI and CT for biopsies containing orthopedic implants was undertaken by Sarve and collaborators [10]. They employed chamfer matching (CM) with only one degree of freedom and therefore, they required the histological specimen to be taken exactly along the implant axis. However, this is a limitation that can hardly be achieved in practice. Their simulated annealing (SA) approach was not initialized with a first guess and used normalized mutual information instead of the binary images to assess similarity. This worked less reliable which was reflected in a high failure rate of 25% [10]. The results of the present study indicate that CM is a solid and fast method to pose a first guess for the bimodal registration of biopsies with dental implants. RMS was chosen as similarity measure to register the initial series of 2D CT slices with HI. The implant edges were aligned successfully also for implants in tilted position in HI. Due to shape similarities of the implants around its axis, however, no clear optima could be found when comparing the 72 RMS scores obtained by rotation of the normal vector around the implant axis. Therefore, additional to the chamfer matching approach described by Borgefors [25], the transformation matrices from CM had to be applied to the segmented bone and the L score could be used to find the optimal orientation of the normal vector. The
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
67
Fig. 8. Regression analysis: the association of 2D HI BIC with 2D CT measures yielded positive regression coefficients for the manual (ˇ = 1.01, p < 0.001) and for the 2 2 = 0.87, a: Rkorr = 0.75). automated (ˇ = 1.7, p < 0.001) measurements (m: Rkorr
subsequent application of SA improved registration results especially for HI specimen taken rather away from the implant axis, such as specimen 6 (Fig. 6). The combined registration with CM and SA using conjunct binary images yielded a highly precise registration with median similarity scores of L = 90%. This L score can hardly be increased closer to hundred percent due to segmentation limitations caused by artifacts and irregular staining of HI specimens or tissue damage on HI production (Fig. 7). The problem is that artifacts can be misclassified and that components from bone may be lost during HI segmentation. However, the registered slices have been compared visually and only minor differences could be observed. For example, in one reconstructed slice there were slight differences in the amount of bone tissue close to the implant (arrow in Fig. 6). Few publications aimed at assessing the agreement of automated BIC measures in CT and manually obtained values from HI. In particular, Vandeweghe et al. found significant correlation between 3D CT BIC and HI BIC [8], but 3D to 2D comparison is problematic because BIC highly depends on the cutting position [27]. Song et al. [28] compared 2D CT and HI BIC on manually registered slices to assess the influence of peri-implant artifacts on bone morphometric analysis. They reported correlation coefficients ranging from 0.29 to 0.43 depending on the distance from the measurement mask to the implant with best correlation at 5 pixel distance. Their results are very different from the present study and the reason may be that manual registration did not respect all degrees of freedom thus giving emphasis on the need for automated approaches. BIC was measured and compared with HI (Fig. 8). The associations were significant both for automated and manual BIC measurements, which additionally confirmed goodness of the registration. While the manually assessed BIC values were almost likewise for both techniques, it was found that the automated BIC values were systematically elevated compared to the manual measurements. This may come from the fact that artifact correction was out of scope for this study and no threshold adoption was performed to segment bone in proximity of the implants. However, the results underline that adoption of thresholds and usage of artifact correction procedures may be suitable to bring manual and automated measurements closer together. The registration
method proposed in this study may be a helpful tool to verify the results. Thus future work should aim to find appropriate artifact correction and threshold adoption strategies to improve bone segmentation in proximity of titanium implants. The resulting corrected BIC values should be validated through direct comparison with the respective HI slide. As soon as both methods yield comparable results, intravital scanners may be used to study osseointegration in longitudinal studies. Hence BIC values obtained in the same animal can be compared and the overall number of animals can be reduced. Furthermore, hard tissue sections may be replaced by paraffin embedded decalcified sections which are much cheaper and better suitable for automated evaluation and immunohistochemistry. Future work should also investigate which extent of mineralization and which trabecular thickness is needed for correct segmentation in proximity of titanium implants. This is essential because in the initial healing phase bone is less mineralized and trabecular thickness is low. Thus differences between HI and CT may be also observed during new bone formation or bone remodeling. However, as the biopsies of this studies were harvested after 20 weeks of healing, mineralization of the bone was already completed. In conclusion, the present work provided a novel reliable and highly accurate method to align corresponding image elements of HI and CT. The visual observations confirmed the high quality of the registration results, which were also corroborated by significant association of BIC measures. However, future research is needed to find suitable thresholds and artifact correction methods for bone segmentation in proximity of titanium implants. This will bring manually obtained BIC values and automatically obtained BIC values closer together and the method provided in this work may be used to verify agreement.
Acknowledgements This work was supported by an unrestricted grant of the Camlog Foundation, Basel, Switzerland. We appreciate the help of Thomas Kaiser, Tina Hagena and Brigitte Hartig who prepared a separate series of histological slices for this study.
68
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compmedimag. 2015.04.005 References [1] Garetto LP, Chen J, Parr JA, Roberts WE. Remodeling dynamics of bone supporting rigidly fixed titanium implants: a histomorphometric comparison in four species including humans. Implant Dent 1995;4:235–43. [2] Liu Y, Enggist L, Kuffer AF, Buser D, Hunziker EB. The influence of BMP-2 and its mode of delivery on the osteoconductivity of implant surfaces during the early phase of osseointegration. Biomaterials 2007;28:2677–86. [3] Kohal RJ, Weng D, Bachle M, Strub JR. Loaded custom-made zirconia and titanium implants show similar osseointegration: an animal experiment. J Periodontol 2004;75:1262–8. [4] Brånemark PI. Osseointegration and its experimental background. J Prosthet Dent 1983;50:399–410. [5] Müller R, Van Campenhout H, Van Damme B, Van der Perre G, Dequeker J, Hildebrand T, et al. Morphometric analysis of human bone biopsies: a quantitative structural comparison of histological sections and micro-computed tomography. Bone 1998;23:59–66. [6] Thomsen JS, Laib A, Koller B, Prohaska S, Mosekilde L, Gowin W. Stereological measures of trabecular bone structure: comparison of 3D micro computed tomography with 2D histological sections in human proximal tibial bone biopsies. J Microsc 2005;218:171–9. [7] Dias DR, Leles CR, Batista AC, Lindh C, Ribeiro-Rotta RF. Agreement between histomorphometry and microcomputed tomography to assess bone microarchitecture of dental implant sites. Clin Implant Dent Relat Res 2013 (epub ahead of print). [8] Vandeweghe S, Coelho PG, Vanhove C, Wennerberg A, Jimbo R. Utilizing micro-computed tomography to evaluate bone structure surrounding dental implants: a comparison with histomorphometry. J Biomed Mater Res B Appl Biomater 2013;101:1259–66. [9] Bernhardt R, Scharnweber D, Müller B, Thurner P, Schliephake H, Wyss P, et al. Comparison of microfocus- and synchrotron X-ray tomography for the analysis of osteointegration around Ti6Al4V implants. Eur Cell Mater 2004;7:42–51, discussion 51. [10] Sarve H, Lindblad J, Johansson CB. Registration of 2D histological images of bone implants with 3D SRmicroCT volumes. In: Bebis G, Boyle RD, Parvin B, Koracin D, Remagnino P, Porikli FM, et al., editors. ISVC. Springer; 2008. p. 1071–80. [11] Ko CC, Douglas WH, DeLong R, Rohrer MD, Swift JQ, Hodges JS, et al. Effects of implant healing time on crestal bone loss of a controlled-load dental implant. J Dent Res 2003;82:585–91. [12] De Smet E, Jaecques SVN, Wevers M, Jansen JA, Jacobs R, Sloten JV, et al. Effect of controlled early implant loading on bone healing and bone mass in guinea pigs, as assessed by micro-CT and histology. Eur J Oral Sci 2006;114:232–42. [13] Xu B, Zhang J, Brewer E, Tu Q, Yu L, Tang J, et al. Osterix enhances BMSCassociated osseointegration of implants. J Dent Res 2009;88:1003–7. [14] Markelj P, Tomaˇzeviˇc D, Likar B, Pernuˇs F. A review of 3D/2D registration methods for image-guided interventions. Med Image Anal 2012;16:642–61. [15] Shalaby A, Farag AA, Ross A, Hockenbury T. 2D–3D registration of human ankle using X-ray and CT images. In: 2012 Cairo Int. Biomed. Eng. Conf. IEEE; 2012. p. 23–6. [16] Fuerst B, Wein W, Müller M, Navab N. Automatic ultrasound–MRI registration for neurosurgery using the 2D and 3D LC2 Metric. Med Image Anal 2014;18:1312–9. [17] Samir C, Kurtek S, Srivastava A, Canis M. Elastic Shape analysis of cylindrical surfaces for 3D/2D registration in endometrial tissue characterization. IEEE Trans Med Imaging 2014;33:1035–43.
[18] Zollei L, Grimson E, Norbash A, Wells W. 2D–3D rigid registration of X-ray fluoroscopy and CT images using mutual information and sparsely sampled histogram estimators. In: Proc. 2001 IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognition. CVPR 2001. IEEE Comput. Soc.; 2001. p. II-696–703. [19] Livyatan H, Yaniv Z, Joskowicz L. Gradient-based 2-D/3-D rigid registration of fluoroscopic X-ray to CT. IEEE Trans Med Imaging 2003;22:1395–406. [20] Varnavas A, Carrell T, Penney G. Fully automated initialisation of 2D–3D image registration. In: 2013 IEEE 10th Int. Symp. Biomed. Imaging. IEEE; 2013. p. 568–71. [21] Birkfellner W, Figl M, Kettenbach J, Hummel J, Homolka P, Schernthaner R, et al. Rigid 2D/3D slice-to-volume registration and its application on fluoroscopic CT images. Med Phys 2007;34:246. [22] Osechinskiy S, Kruggel F. Slice-to-volume nonrigid registration of histological sections to MR images of the human brain. Anat Res Int 2011;2011. [23] Thibault GP, Stack R, Sesterhenn IA, Potter K, Meyer CR. Registration methods for histological slides and ex vivo MRI of prostate. In: 2007 IEEE Nucl. Sci. Symp. Conf. Rec. IEEE; 2007. p. 3102–4. [24] Schwarz F, Mihatovic I, Golubovich V, Schär A, Sager M, Becker J. Impact of abutment microstructure and insertion depth on crestal bone changes at nonsubmerged titanium implants with platform switch. Clin Oral Implants Res 2014 (epub ahead of print). [25] Borgefors G. Hierarchical chamfer matching: a parametric edge matching algorithm. IEEE Trans Pattern Anal Mach Intell 1988;10:849–65. [26] Development Core Team R. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2011. [27] Bernhardt R, Kuhlisch E, Schulz MC, Eckelt U, Stadlinger B. Comparison of bone-implant contact and bone-implant volume between 2D-histological sections and 3D-SRmicroCT slices. Eur Cell Mater 2012;23:237–47, discussion 247–248. [28] Song JW, Cha JY, Bechtold TE, Park YC. Influence of peri-implant artifacts on bone morphometric analysis with micro-computed tomography. Int J Oral Maxillofac Implants 2013;28:519–25. Dr. Kathrin Becker studied dentistry during 2005–2013 and applied com¨ (Germany). Since puter science during 2007–2014 in Greifswald and Gottingen November 2012, she is an elected member of the DGI committee Nexte Generation. From February 2014 she has been a Research assistant (DFG project) at Heinrich-Heine University, Düsseldorf. Her research interests are on microCT, image processing, machine learning methods, implant surfaces and osseointegration. Dr. Martin Stauber studied physics during 1995–2000, had done P.G. course in Medical Physics in 2000–2002, did PhD in 2005 and pursued Postdoc at the Institute for Biomechanics, ETH Zürich in 2007. He worked as Lecturer during 2007–2012 as well as a Senior Scientist at Scanco Medical AG from 2012. His research interests are on microCT, image processing, bone biology, bone microstructure, and bone mechanics. Prof. Dr. Frank Schwarz studied dentistry at University of Saarland, Homburg, Germany during 1993–1998 and continued with P.G in 1999–2000 in the same institution. He worked as Research Associate between 2000 and 2002, as an Assistant Professor between 2002 and 2006 and from 2006 as an Associate Professor. Besides his editorial works, he is the Vice President German Association of Oral Implantology (DGI). His research interests are on guided bone regeneration, growth factors, implant surfaces, and peri-implant diseases. Prof. Dr. Tim Beißbarth studied Biology and Computer Science during 1992–1998 and completed PhD in 2001. His academic life started as Postdoc in 2001–2002 at Berlin, in 2002–2005 in Melbourne and as GL at Heidelberg for various institutions. Currently, he is an Associate Professor in University Medical Center Göttingen. His research interests are on bioinformatics, systems biology, machine learning methods, analysis of high-throughput data, oncology, ontologies and biological signaling pathways.
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Automated 3D–2D registration of X-ray microcomputed tomography with histological sections for dental implants in bone using chamfer matching and simulated annealing Kathrin Becker a,b,∗ , Martin Stauber c , Frank Schwarz b , Tim Beißbarth a a
Department of Medical Statistics, Biostatistics Group, University Medical Center, Georg-August University, Humboldt Allee 32, 37073 Göttingen, Germany Department of Oral Surgery, Westdeutsche Kieferklinik, Heinrich-Heine University, Moorenstr. 5, 40225 Düsseldorf, Germany c Scanco Medical AG, Fabrikweg 2, 8306 Brüttisellen, Switzerland b
a r t i c l e
i n f o
Article history: Received 15 November 2014 Received in revised form 1 March 2015 Accepted 17 April 2015 JEL classification: I10 C65 Keywords: Image registration Micro-computed tomography Histology Osseointegration Dental implants Histomorphometry
a b s t r a c t We propose a novel 3D–2D registration approach for micro-computed tomography (CT) and histology (HI), constructed for dental implant biopsies, that finds the position and normal vector of the oblique slice from CT that corresponds to HI. During image pre-processing, the implants and the bone tissue are segmented using a combination of thresholding, morphological filters and component labeling. After this, chamfer matching is employed to register the implant edges and fine registration of the bone tissues is achieved using simulated annealing. The method was tested on n = 10 biopsies, obtained at 20 weeks after non-submerged healing in the canine mandible. The specimens were scanned with CT 100 and processed for hard tissue sectioning. After registration, we assessed the agreement of bone to implant contact (BIC) using automated and manual measurements. Statistical analysis was conducted to test the agreement of the BIC measurements in the registered samples. Registration was successful for all specimens and agreement of the respective binary images was high (median: 0.90, 1.–3. Qu.: 0.89–0.91). Direct comparison of BIC yielded that automated (median 0.82, 1.–3. Qu.: 0.75–0.85) and manual (median 0.61, 1.–3. Qu.: 0.52–0.67) measures from CT were significant positively correlated with HI (median 0.65, 1.–3. Qu.: 0.59–0.72) between CT and HI groups (manual: R2 = 0.87, automated: R2 = 0.75, p < 0.001). The results show that this method yields promising results and that CT may become a valid alternative to assess osseointegration in three dimensions. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Biological processes associated with the osseointegration of endosseous dental implants have been extensively investigated in numerous preclinical and clinical studies [1–3]. The level of osseointegration is typically evaluated by histological assessment of bone to implant contact (BIC) areas in undecalcified hard tissue sections [4]. One major drawback of nondecalcified sectioning, however, is a significant loss of biological
∗ Corresponding author at: Department of Orthodontics, Westdeutsche Kieferklinik, Heinrich-Heine University, Moorenstr. 5, 40225 Düsseldorf, Germany. Tel.: +49 211 8118140; fax: +49 211 8116550. E-mail addresses: [email protected] (K. Becker), [email protected] (M. Stauber), [email protected] (F. Schwarz), [email protected] (T. Beißbarth). http://dx.doi.org/10.1016/j.compmedimag.2015.04.005 0895-6111/© 2015 Elsevier Ltd. All rights reserved.
tissues during the histological (HI) processing of serial sections, and therefore, only a maximum of 2–3 sections are available which may not reflect the entire healing processes at the investigated site. Another drawback is the amount of time needed for the processing of hard tissue sections, delaying data acquisition for several months after completion of the experimental procedures. Moreover, staining of specimens is a highly sensitive process and commonly associated with marked variances in color intensity. Furthermore, evaluation is limited to 2D, so different orientations can not be investigated, and killing of the animal is needed to harvest biopsies for histology thus for each time point different animals are needed. In recent years, micro computed tomography (CT) has become a promising tool for quantitative three-dimensional (3D) analyses of new bone formation following implant placement and various types of bone augmentation procedures. This non-destructive and non-invasive volumetric analysis offers many advantages including
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
a fast and automated characterization of structural features such as bone volume density [5,6] or the assessment of BIC values [7,8] in different planes and even in 3D. Furthermore, with the advent of intravital CT scanners, different time points can be studied within the same animal. Despite these advantages, assessment of BIC values in CT scans may be impeded by beam hardening artifacts, partial volume effects and scattering [9,10] which occur in the proximity of dental titanium implants. While partial volume effects can be eliminated by performing BIC measurements few voxels away from the implant, beam hardening and scatter artifacts may distort BIC measurements. Due to the lacking proof of comparability, CT analyses are usually complemented with conventional HI at the moment [11–13]. To verify if CT and HI yield the same results, and to validate novel BIC measurements in the future, 2D corresponding slices have to be compared with the gold standard HI. For this a 3D–2D registration is required that aligns corresponding image elements from CT and HI. Few approaches for 3D–2D registration are described in literature, mostly reporting on the registration of CT or MRI volumes with projected 2D X-ray [14,15], ultrasound [16,17] or fluoroscopic images [18–20]. Registration of single 2D slices obtained during interventional procedures with CT [21] or HI section with CT did not gain much attention so far since the registration relies on manyto-one correspondence of the underlying 3D-to-2D mapping as well as a statistical consistency problem owing to the reliance of a similarity measure on small sampling size of an input 2D image [22] and big differences in information content [23]. Manual approaches for 3D–2D registration are hard to realize because HI is usually not located parallel to a reference plane within CT, thus, a seven dimensional optimization problem needs to be solved (rotation, translation and scaling). Sarve et al. were the first to report on automated computerized registration for synchroton radiation microcomputed tomography and histology using the model of osseointegrated orthopedic titanium (TI) implants [10]. They used normalized mutual information to assess similarity and assumed HI being located directly at the implant axis. They also tried to employ simulated annealing without a prior first guess and fast cooling so a failure rate of 25% was reported by the authors [10]. The aim of the present study was to develop a robust automated computerized registration method for HI images with CT volumes. The specific goals were (i) to establish an algorithm to automatically register 2D HI images with a 3D CT volume, and (ii) to assess the agreement of BIC values in the registered images.
2. Acquisition of data 2.1. Tissue samples The specimens (n = 10) used for the present analysis were randomly selected from tissue biopsies obtained in a recently published experimental study [24]. In brief, each biopsy contained one commercially available titanium implant (Ø 3.8 mm, length: 9 mm, Conelog® Screw-Line Implant, Camlog Biotechnologies AG, Basel, Switzerland) connected with a diameter-reduced wide body Ti (Ti6Al4V) healing abutment (Ø 3.8 mm, height: 4 mm, Camlog). These biopsies were retrieved at 20 weeks after implant placement and non-submerged healing in the mandible of foxhounds. The study protocol was approved by the appropriate local authority (Landesamt fr Natur und Verbraucherschutz, Recklinghausen, Germany). The biopsies were dehydrated using ascending grades of alcohol and xylene, infiltrated and embedded in methylmethacrylate (Technovit 9100 NEU; Heraeus Kulzer, Wehrheim, Germany),
63
which is the standardized procedure for non-decalcified sectioning. This procedure takes 20 h, during which negative influence of polymerization heat was avoided by storing the specimen at cold atmosphere (4 ◦ C). 2.2. MicroCT acquisition After embedding in methylmethacrylate, the biopsies were scanned with a CT 100 (Scanco Medical AG, Brttisellen, Switzerland). The X-ray tube was operated at an energy of 90 kVp and an intensity of 88 A. The integration time was chosen to 1200 ms and the frame averaging was set at 2×. The data was reconstructed into 3D volumes with an isotropic nominal resolution of 8.6 m. 2.3. Histological preparation and image acquisition The technique has been described in detail [24]. Briefly, after CT analysis, histological preparation was performed. Each implant site was cut in the vestobulo-oral direction along with the long axis of the implant using a diamond band saw (Exact(r); Apparatebau, Nordersted, Germany). Serial sections of approximately 300 m in thickness were prepared, glued with acrylic cement (Technovit 7210 VLC; Hareus Kulzer), to silanized glass slides (Super Frost; Menzel GmbH, Braunschweig, Germany) and ground to a final thickness of approximately 40 m. Subsequently, from each tissue block, one specimen was randomly selected for the present analysis and stained with paragon. For image acquisition a color CCD camera (Color View III, Olympus, Hamburg, Germany) mounted on a binocular light microscope (Olympus BX50, Olympus, Hamburg, Germany) was used. Digital images for registration (original magnification 100×) were obtained using the Cell D program (Cell D® , Soft Imaging System, Münster, Germany). Additionally, the specimens were scanned with Aperio ScanScope Slide Scanner (Aperio Technologies Inc., CA, USA) for BIC assessment (original magnification 200×) and the region of interest was exported as snapshot using Aperio ImageScope v11.1.2.752 (Aperio Technologies). 3. Image registration 3.1. Image preprocessing and feature extraction The image preprocessing and registration was performed with Matlab (Matlab R2013a 64-bit, The Mathworks Inc., Massachusetts, USA). During preprocessing, both HI and CT were each transformed into two binary images showing either the dental implant or the bone (Fig. 1). To segment the implant from HI, thresholding and component labeling was applied. The edges were detected using Canny edge filter. For segmentation of the bone tissue from HI, a combination of thresholding, component labeling and morphological filters (erosion, dilation) with a 3 × 3 disk as structuring element was used. On the CT slices, the thresholds were chosen manually such that either the bone tissue or dental implant was segmented. 3.2. Adaption of image resolution (scaling) The element size of histology was obtained from a bar of fixed length (5 mm) located within the histological images. The CT element size was read from the header of the CT-files given in Scanco Medical AIM-format. Adaption was achieved by scaling HI images with the ratio of both element sizes, i.e. sizehi /sizect .
64
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Fig. 1. Image preprocessing procedure: to segment the bone tissue from HI, a mask to cover the bone was computed using thresholding and morphological filters. With component labeling, the marrow spaces could be cleared. Segmentation of the implant from HI was conducted by finding the largest component. At the CT images segmentation was done using thresholding.
the approximated Euclidean distance to the nearest edge, the alignment score was given by the root mean squared average (RMS), which averages the intensities xi to which the ei were mapped:
n 11 RMS = xi n
3
(1)
i=1
The division by 3 is required since the 3–4 Chamfer distance was used instead of the Euclidean distance. The optimal match was found when no transformations were possible that would reduce RMS. Convergence of CM was shown by Borgefors [25]. Due to proximate rotation symmetry of the implant, RMS from CM with slices located around the implant axis was expected to be likewise. Thus the segmented bone from HI was mapped to the segmented bone from CT slice using the transformation matrix from CM. Then, similarity L was measured between the binary images by averaging the L1 norm and subtraction from 1 (xi , yi denote the values of aligned pixels): Fig. 2. Chamfer matching principle: the distance transform (DT) of 72 slices, that had been extracted from the 3D volume before and were located around the implant axis, was computed. The output was the optimal alignment of the corresponding 2D edges.
3.3. First guess registration using chamfer matching and principal component analysis Chamfer matching (CM) was used to find a first guess for the start of the fine registration and was implemented as described by Borgefors [25] for the 2D case (Fig. 2). Two pyramid levels were used and initialization was achieved by finding the implant axis from CT, which was approximated with the first axis (PA) from principal component analysis (PCA): to conduct PCA, all edge vectors from HI and CT were centered and stored as row vectors vi in a matrix. Then, its covariance matrix and eigenvalue decomposition were computed. The eigenvector with largest eigenvalue yielded PA. Next, a vector orthogonal to the implant axis (third principal axis, i.e. eigenvector with smallest eigenvalue) was rotated in 5 degree steps around PA and a series of 72 slices was extracted from the CT volume. From each slice, the 3–4 distance transform (DT) was computed [25] and the implant edges ei from HI were moved along DT until an optimal match was found. Since each intensity of DT reflects
1 xi − yi 1 n n
L =1−
(2)
i=1
Hence, RMS was used to find the transformation of HI and a slice extracted from CT. L was used to compare similarity of different slices extracted from CT with HI. These slices were obtained by rotating a vector orthogonal to the implant axis in 5 degree steps. At the position with highest L score, additional 1 degree rotations of the normal were performed (10◦ to each direction). Missing intensity values of oblique slices were interpolated using trilinear interpolation. The parameters revealing the optimal L score finally yielded the first guess registration of the implants using CM. 3.4. Fine registration of the segmented bone using simulated annealing Simulated annealing (SA) is a robust optimization method incorporated from thermodynamics. It was implemented with 6 degrees of freedom (rotation, translation) to find a slice from CT that maximizes L when being aligned with HI. Input were the segmented bone binary images from CT and HI as well as the so far optimal position and normal vector from Chamfer Matching. 3D–2D
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Fig. 3. Manual BIC measurements: a segmented line tool was used to measure the lengths of the implant outline (left) and the implant fractions being in direct contact with bone (right, green arrow) below the implant shoulders. The sum of the fractions divided by the length of the implant outline was defined as BIC measure. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
65
Fig. 4. Automated BIC measurements: after thresholding (right image), a mask was located 3 pixel distant from the implant outline. The number of pixels classified as bone (red pixels) divided by the entire number of pixels within the mask (red or green pixels) yielded the automated 2D BIC measure. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
matching was achieved by refining the position and normal vector of the CT slice, while HI remained fixed. The algorithm uses a parameter temperature T that was initialized with 300 and cooled down to 30. The step width for transformations was set to s = T/100. The Boltzmann constant k was set to 0.01. The energy states E1 , E2 were set to the alignment scores L of the present and next slice to be extracted from the volume. The next transformation was chosen randomly and excepted if the probability p of a system to change energy from E1 to E2 was smaller than a random number. p=
⎧
⎨ exp − E1 − E2 if E2 ≤ E1 ⎩
kT
1
if E2 > E1
The temperature was cooled down by 1% after 250 iterations and parameters yielding greatest L finally provided the optimal alignment. The code of the implementation can be found at: http:// kathrin-becker.eu/code/Registration.zip
Fig. 5. Threshold intervals were selected such that the implant as well as bone tissue could be segmented clearly from CT. However, in proximity to the implant (red arrow, right image) artifacts were sometimes mis-classified as bone (green arrow, left image) and thus increased automated BIC measures (colored mask, left image). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
3.5. Quantification of registration results and BIC measurements For quantification of the registration results, the L scores from registration were used. In addition, the registration results were visually examined and to evaluate agreement of BIC between the two imaging modalities, BIC was assessed in the region below the implant shoulder both automated and manually. This measurement refers to the relative fraction of bone being in direct contact with the implant. Manual measurements were performed using ImageJ (National Institutes of Health, USA) by two experienced investigators (histology: KB, CT: MS). The length of the implant contour and the lengths of the implant outline in direct contact with the bone were measured in the region below the implant shoulders using the Segmented Line measurement tool (Fig. 3). Prior to the start of the analyses, a calibration procedure was initiated for the investigators and revealed that repeated measurements were similar at >95% level. Automated measurements to obtain CT BIC were implemented and performed using Matlab R2013a. The automated BIC was assessed within a mask of 2 pixel width located below the implant shoulder and 3 pixel distant from the implant surface (Fig. 4). The mask was computed by applying image dilation and erosion filters
on the segmented implant. The position of the implant shoulders was specified manually using Matlab’s impixel feature. The BIC measure was given by the relative fraction of pixels classified as bone. To segment the bone, the same procedure as described in Section 3.1 was used. 3.6. Statistical analysis The statistical analyses were performed using the R statistics software version 2.12.2 [26]. Test for linear association between 2D HI BIC and 2D CT BIC was performed by calculation of linear models. Results were considered statistically significant at p < 0.05. 4. Results 4.1. Image preprocessing Segmentation of the implants from CT was successful for all specimens investigated. Segmentation of bone tissue from CT was sometimes impeded by artifacts, which were classified as bone in close proximity to the implant surface (Fig. 5).
66
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Table 1 The registration results are presented column-wise. Similarity scores L: chamfer matching (CM), simulated annealing (SA). BIC measurements: histology (HI), manual CT measurements (m), automated CT BIC measurements (a). The last rows shows the respective mean values ± standard deviation, median and quartiles. Specimen
L (CM)
L (SA)
HI BIC
CT BIC(m)
CT BIC(a)
1 2 3 4 5 6 7 8 9 10
0.86 0.83 0.83 0.88 0.83 0.87 0.83 0.80 0.87 0.82
0.91 0.90 0.93 0.93 0.90 0.91 0.85 0.87 0.89 0.92
0.60 0.77 0.73 0.80 0.46 0.64 0.65 0.59 0.22 0.72
0.61 0.79 0.68 0.75 0.52 0.61 0.55 0.50 0.26 0.62
0.86 0.82 0.86 0.91 0.75 0.89 0.77 0.74 0.61 0.81
Mean ± SD Median 1.–3. Qu.
0.85 ± 0.03 0.83 0.83–0.87
0.9 ± 0.03 0.91 0.9–0.92
0.62 ± 0.17 0.65 0.59–0.73
0.59 ± 0.15 0.61 0.52–0.67
0.8 ± 0.09 0.82 0.75–0.85
Fig. 7. The registration of specimen 7 was reliable even though parts of the bone tissue were damaged on production and despite visible color variance of the staining agent. However, this was reflected in the lowest L score of 85 %. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
4.3. Assessment of bone to implant contact Significant correlation of BIC for both, manual (Pearson = 0.94, p < 0.001, CI: 0.80–1) and automated measurements were found (Pearson = 0.88, p < 0.001, CI: 0.65–1). Regression analysis yielded a significant association between HI and CT groups (manual mea2 surements: regression coefficient ˇ = 1.01, Rkorr = 0.87, p < 0.001; 2 = automated measurements: regression coefficient ˇ = 1.7, Rkorr 0.75, p < 0.001). BIC values were higher for the automated measurements (median: 0.82; 1.–3. Qu.: 0.75–0.85) compared to the manual measurements (median 0.61; 1.–3. Qu.: 0.52*-0.67) and HI (median 0.65; 1.–3. Qu.: 0.59–0.72). Artifacts were present mainly in the apical portion of the CT slices. Fig. 6. Registration results of specimen 6: this specimen was not sawed directly at its center thus chamfer matching could only provide a rough approximation. Interestingly, the marrow space at the arrow is hardly visible in CT.
The implants could be segmented clearly from HI without adjustment of threshold levels. The bone tissue segmentation from HI required minor adjustments of threshold levels due to low staining intensities in three specimens. As small components were removed during segmentation, very small bone elements were missing after segmentation. It was found that this segmentation procedure was very sensitive to the staining process and may need to be adjusted for different staining procedures.
4.2. Image registration Computation of the implant axes and subsequent application of chamfer matching provided reliable first guesses when using similarity of the segmented bone L for inter-slice comparison. For the 2D matching of the segmented bone, RMS was in the range of 0.47–4.10. When comparing the RMS values obtained by rotating the normal around the implant axes, no clear global optimum could be found for RMS scores among different slices due to shape similarities of the implants around their axes. The transformation matrix from CM was used to match also the bone tissue which, in contrast to RMS, yielded clear optima and thus a reliant first guess for SA. The similarity scores L were on median 0.83 with 1st to 3rd quartile reaching from 0.83 to 0.87 (Table 1). Simulated annealing provided successful fine registration for all specimen and increased the similarity scores L in all cases (median 0.9, 1.–3. Qu.: 0.89–0.91) (Table 1). Visual examination corroborated the high agreement of marrow spaces and bone crest (Fig. 6).
5. Discussion In this study, we presented a novel method for automated 3D–2D registration of X-ray microcomputed tomography with histological sections. This method was developed for the direct comparison of these two techniques, and to evaluate the agreement of the bone to implant contact (BIC). To prove that the correct position and normal vector for the corresponding CT slice were found, the results were quantified by means of similarity scores and by visual examination. To the best of the authors’ knowledge, the only study dealing with registration of HI and CT for biopsies containing orthopedic implants was undertaken by Sarve and collaborators [10]. They employed chamfer matching (CM) with only one degree of freedom and therefore, they required the histological specimen to be taken exactly along the implant axis. However, this is a limitation that can hardly be achieved in practice. Their simulated annealing (SA) approach was not initialized with a first guess and used normalized mutual information instead of the binary images to assess similarity. This worked less reliable which was reflected in a high failure rate of 25% [10]. The results of the present study indicate that CM is a solid and fast method to pose a first guess for the bimodal registration of biopsies with dental implants. RMS was chosen as similarity measure to register the initial series of 2D CT slices with HI. The implant edges were aligned successfully also for implants in tilted position in HI. Due to shape similarities of the implants around its axis, however, no clear optima could be found when comparing the 72 RMS scores obtained by rotation of the normal vector around the implant axis. Therefore, additional to the chamfer matching approach described by Borgefors [25], the transformation matrices from CM had to be applied to the segmented bone and the L score could be used to find the optimal orientation of the normal vector. The
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
67
Fig. 8. Regression analysis: the association of 2D HI BIC with 2D CT measures yielded positive regression coefficients for the manual (ˇ = 1.01, p < 0.001) and for the 2 2 = 0.87, a: Rkorr = 0.75). automated (ˇ = 1.7, p < 0.001) measurements (m: Rkorr
subsequent application of SA improved registration results especially for HI specimen taken rather away from the implant axis, such as specimen 6 (Fig. 6). The combined registration with CM and SA using conjunct binary images yielded a highly precise registration with median similarity scores of L = 90%. This L score can hardly be increased closer to hundred percent due to segmentation limitations caused by artifacts and irregular staining of HI specimens or tissue damage on HI production (Fig. 7). The problem is that artifacts can be misclassified and that components from bone may be lost during HI segmentation. However, the registered slices have been compared visually and only minor differences could be observed. For example, in one reconstructed slice there were slight differences in the amount of bone tissue close to the implant (arrow in Fig. 6). Few publications aimed at assessing the agreement of automated BIC measures in CT and manually obtained values from HI. In particular, Vandeweghe et al. found significant correlation between 3D CT BIC and HI BIC [8], but 3D to 2D comparison is problematic because BIC highly depends on the cutting position [27]. Song et al. [28] compared 2D CT and HI BIC on manually registered slices to assess the influence of peri-implant artifacts on bone morphometric analysis. They reported correlation coefficients ranging from 0.29 to 0.43 depending on the distance from the measurement mask to the implant with best correlation at 5 pixel distance. Their results are very different from the present study and the reason may be that manual registration did not respect all degrees of freedom thus giving emphasis on the need for automated approaches. BIC was measured and compared with HI (Fig. 8). The associations were significant both for automated and manual BIC measurements, which additionally confirmed goodness of the registration. While the manually assessed BIC values were almost likewise for both techniques, it was found that the automated BIC values were systematically elevated compared to the manual measurements. This may come from the fact that artifact correction was out of scope for this study and no threshold adoption was performed to segment bone in proximity of the implants. However, the results underline that adoption of thresholds and usage of artifact correction procedures may be suitable to bring manual and automated measurements closer together. The registration
method proposed in this study may be a helpful tool to verify the results. Thus future work should aim to find appropriate artifact correction and threshold adoption strategies to improve bone segmentation in proximity of titanium implants. The resulting corrected BIC values should be validated through direct comparison with the respective HI slide. As soon as both methods yield comparable results, intravital scanners may be used to study osseointegration in longitudinal studies. Hence BIC values obtained in the same animal can be compared and the overall number of animals can be reduced. Furthermore, hard tissue sections may be replaced by paraffin embedded decalcified sections which are much cheaper and better suitable for automated evaluation and immunohistochemistry. Future work should also investigate which extent of mineralization and which trabecular thickness is needed for correct segmentation in proximity of titanium implants. This is essential because in the initial healing phase bone is less mineralized and trabecular thickness is low. Thus differences between HI and CT may be also observed during new bone formation or bone remodeling. However, as the biopsies of this studies were harvested after 20 weeks of healing, mineralization of the bone was already completed. In conclusion, the present work provided a novel reliable and highly accurate method to align corresponding image elements of HI and CT. The visual observations confirmed the high quality of the registration results, which were also corroborated by significant association of BIC measures. However, future research is needed to find suitable thresholds and artifact correction methods for bone segmentation in proximity of titanium implants. This will bring manually obtained BIC values and automatically obtained BIC values closer together and the method provided in this work may be used to verify agreement.
Acknowledgements This work was supported by an unrestricted grant of the Camlog Foundation, Basel, Switzerland. We appreciate the help of Thomas Kaiser, Tina Hagena and Brigitte Hartig who prepared a separate series of histological slices for this study.
68
K. Becker et al. / Computerized Medical Imaging and Graphics 44 (2015) 62–68
Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compmedimag. 2015.04.005 References [1] Garetto LP, Chen J, Parr JA, Roberts WE. Remodeling dynamics of bone supporting rigidly fixed titanium implants: a histomorphometric comparison in four species including humans. Implant Dent 1995;4:235–43. [2] Liu Y, Enggist L, Kuffer AF, Buser D, Hunziker EB. The influence of BMP-2 and its mode of delivery on the osteoconductivity of implant surfaces during the early phase of osseointegration. Biomaterials 2007;28:2677–86. [3] Kohal RJ, Weng D, Bachle M, Strub JR. Loaded custom-made zirconia and titanium implants show similar osseointegration: an animal experiment. J Periodontol 2004;75:1262–8. [4] Brånemark PI. Osseointegration and its experimental background. J Prosthet Dent 1983;50:399–410. [5] Müller R, Van Campenhout H, Van Damme B, Van der Perre G, Dequeker J, Hildebrand T, et al. Morphometric analysis of human bone biopsies: a quantitative structural comparison of histological sections and micro-computed tomography. Bone 1998;23:59–66. [6] Thomsen JS, Laib A, Koller B, Prohaska S, Mosekilde L, Gowin W. Stereological measures of trabecular bone structure: comparison of 3D micro computed tomography with 2D histological sections in human proximal tibial bone biopsies. J Microsc 2005;218:171–9. [7] Dias DR, Leles CR, Batista AC, Lindh C, Ribeiro-Rotta RF. Agreement between histomorphometry and microcomputed tomography to assess bone microarchitecture of dental implant sites. Clin Implant Dent Relat Res 2013 (epub ahead of print). [8] Vandeweghe S, Coelho PG, Vanhove C, Wennerberg A, Jimbo R. Utilizing micro-computed tomography to evaluate bone structure surrounding dental implants: a comparison with histomorphometry. J Biomed Mater Res B Appl Biomater 2013;101:1259–66. [9] Bernhardt R, Scharnweber D, Müller B, Thurner P, Schliephake H, Wyss P, et al. Comparison of microfocus- and synchrotron X-ray tomography for the analysis of osteointegration around Ti6Al4V implants. Eur Cell Mater 2004;7:42–51, discussion 51. [10] Sarve H, Lindblad J, Johansson CB. Registration of 2D histological images of bone implants with 3D SRmicroCT volumes. In: Bebis G, Boyle RD, Parvin B, Koracin D, Remagnino P, Porikli FM, et al., editors. ISVC. Springer; 2008. p. 1071–80. [11] Ko CC, Douglas WH, DeLong R, Rohrer MD, Swift JQ, Hodges JS, et al. Effects of implant healing time on crestal bone loss of a controlled-load dental implant. J Dent Res 2003;82:585–91. [12] De Smet E, Jaecques SVN, Wevers M, Jansen JA, Jacobs R, Sloten JV, et al. Effect of controlled early implant loading on bone healing and bone mass in guinea pigs, as assessed by micro-CT and histology. Eur J Oral Sci 2006;114:232–42. [13] Xu B, Zhang J, Brewer E, Tu Q, Yu L, Tang J, et al. Osterix enhances BMSCassociated osseointegration of implants. J Dent Res 2009;88:1003–7. [14] Markelj P, Tomaˇzeviˇc D, Likar B, Pernuˇs F. A review of 3D/2D registration methods for image-guided interventions. Med Image Anal 2012;16:642–61. [15] Shalaby A, Farag AA, Ross A, Hockenbury T. 2D–3D registration of human ankle using X-ray and CT images. In: 2012 Cairo Int. Biomed. Eng. Conf. IEEE; 2012. p. 23–6. [16] Fuerst B, Wein W, Müller M, Navab N. Automatic ultrasound–MRI registration for neurosurgery using the 2D and 3D LC2 Metric. Med Image Anal 2014;18:1312–9. [17] Samir C, Kurtek S, Srivastava A, Canis M. Elastic Shape analysis of cylindrical surfaces for 3D/2D registration in endometrial tissue characterization. IEEE Trans Med Imaging 2014;33:1035–43.
[18] Zollei L, Grimson E, Norbash A, Wells W. 2D–3D rigid registration of X-ray fluoroscopy and CT images using mutual information and sparsely sampled histogram estimators. In: Proc. 2001 IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognition. CVPR 2001. IEEE Comput. Soc.; 2001. p. II-696–703. [19] Livyatan H, Yaniv Z, Joskowicz L. Gradient-based 2-D/3-D rigid registration of fluoroscopic X-ray to CT. IEEE Trans Med Imaging 2003;22:1395–406. [20] Varnavas A, Carrell T, Penney G. Fully automated initialisation of 2D–3D image registration. In: 2013 IEEE 10th Int. Symp. Biomed. Imaging. IEEE; 2013. p. 568–71. [21] Birkfellner W, Figl M, Kettenbach J, Hummel J, Homolka P, Schernthaner R, et al. Rigid 2D/3D slice-to-volume registration and its application on fluoroscopic CT images. Med Phys 2007;34:246. [22] Osechinskiy S, Kruggel F. Slice-to-volume nonrigid registration of histological sections to MR images of the human brain. Anat Res Int 2011;2011. [23] Thibault GP, Stack R, Sesterhenn IA, Potter K, Meyer CR. Registration methods for histological slides and ex vivo MRI of prostate. In: 2007 IEEE Nucl. Sci. Symp. Conf. Rec. IEEE; 2007. p. 3102–4. [24] Schwarz F, Mihatovic I, Golubovich V, Schär A, Sager M, Becker J. Impact of abutment microstructure and insertion depth on crestal bone changes at nonsubmerged titanium implants with platform switch. Clin Oral Implants Res 2014 (epub ahead of print). [25] Borgefors G. Hierarchical chamfer matching: a parametric edge matching algorithm. IEEE Trans Pattern Anal Mach Intell 1988;10:849–65. [26] Development Core Team R. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2011. [27] Bernhardt R, Kuhlisch E, Schulz MC, Eckelt U, Stadlinger B. Comparison of bone-implant contact and bone-implant volume between 2D-histological sections and 3D-SRmicroCT slices. Eur Cell Mater 2012;23:237–47, discussion 247–248. [28] Song JW, Cha JY, Bechtold TE, Park YC. Influence of peri-implant artifacts on bone morphometric analysis with micro-computed tomography. Int J Oral Maxillofac Implants 2013;28:519–25. Dr. Kathrin Becker studied dentistry during 2005–2013 and applied com¨ (Germany). Since puter science during 2007–2014 in Greifswald and Gottingen November 2012, she is an elected member of the DGI committee Nexte Generation. From February 2014 she has been a Research assistant (DFG project) at Heinrich-Heine University, Düsseldorf. Her research interests are on microCT, image processing, machine learning methods, implant surfaces and osseointegration. Dr. Martin Stauber studied physics during 1995–2000, had done P.G. course in Medical Physics in 2000–2002, did PhD in 2005 and pursued Postdoc at the Institute for Biomechanics, ETH Zürich in 2007. He worked as Lecturer during 2007–2012 as well as a Senior Scientist at Scanco Medical AG from 2012. His research interests are on microCT, image processing, bone biology, bone microstructure, and bone mechanics. Prof. Dr. Frank Schwarz studied dentistry at University of Saarland, Homburg, Germany during 1993–1998 and continued with P.G in 1999–2000 in the same institution. He worked as Research Associate between 2000 and 2002, as an Assistant Professor between 2002 and 2006 and from 2006 as an Associate Professor. Besides his editorial works, he is the Vice President German Association of Oral Implantology (DGI). His research interests are on guided bone regeneration, growth factors, implant surfaces, and peri-implant diseases. Prof. Dr. Tim Beißbarth studied Biology and Computer Science during 1992–1998 and completed PhD in 2001. His academic life started as Postdoc in 2001–2002 at Berlin, in 2002–2005 in Melbourne and as GL at Heidelberg for various institutions. Currently, he is an Associate Professor in University Medical Center Göttingen. His research interests are on bioinformatics, systems biology, machine learning methods, analysis of high-throughput data, oncology, ontologies and biological signaling pathways.
