- Email: [email protected]
Fast and accurate automated measurements in digitized stereophotogrammetric radiographs
Journal of Biomechanics 31 (1998) 491 — 498
Technical Note
Fast and accurate automated measurements in digitized stereophotogrammetric radiographs Henri A. Vrooman!,*, Edward R. Valstar", Gert-Jan Brand!, Dennis R. Admiraal!, Piet M. Rozing", Johan H.C. Reiber! ! Laboratory for Clinical and Experimental Image Processing, Division of Image processing, Department of Radiology, Leiden University Medical Center, P.O. Box 9600, NL-2300 RC Leiden, The Netherlands " Department of Orthopaedics, Leiden University Medical Center, The Netherlands Received in final form 16 January 1998
Abstract Until recently, Roentgen Stereophotogrammetric Analysis (RSA) required the manual definition of all markers using a highresolution measurement table. To automate this tedious and time-consuming process and to eliminate observer variabilities, an analytical software package has been developed and validated for the detection, identification, and matching of markers in RSA radiographs. The digital analysis procedure consisted of the following steps: (1) the detection of markers using a variant of the Hough circle-finder technique; (2) the identification and labeling of the detected markers; (3) the reconstruction of the three-dimensional position of the bone markers and the prosthetic markers; and (4) the computation of micromotion. To assess the influence of film digitization, the measurements obtained from nine phantom radiographs using two different film scanners were compared with the results obtained by manual processing. All markers in the phantom radiographs were automatically detected and correctly labeled. The best results were obtained with a Vidar VXR-12 CCD scanner, for which the measurement errors were comparable to the errors associated with the manual approach. To assess the in vivo reproducibility, 30 patient radiographs were analyzed twice with the manual as well as with the automated procedure. Approximately, 85% of all calibration markers and bone markers were automatically detected and correctly matched. The calibration errors and the rigid-body errors show that the accuracy of the automated procedure is comparable to the accuracy of the manual procedure. The rigid-body errors had comparable mean values for both techniques: 0.05 mm for the tibia and 0.06 mm for the prosthesis. The reproducibility of the automated procedure showed to be slightly better than that of the manual procedure. The maximum errors in the computed translation and rotation of the tibial component were 0.11 mm and 0.24°, compared to 0.13 mm and 0.27° for the manual RSA procedure. The total processing time is less than 10 min per radiograph, including interactive corrections, compared to approximately 1 h for the manual approach. In conclusion, a new and widely applicable, computer-assisted technique has become available to detect, identify, and match markers in RSA radiographs and to assess the micromotion of endoprostheses. This new technique will be used in our clinic for our hip, knee, and elbow studies. ( 1998 Elsevier Science Ltd. All rights reserved. Keywords: Endoprostheses; Micromotion; Roentgen stereophotogrammetric analysis; Film digitization; Digital image processing
1. Introduction At the Department of Orthopaedics, Leiden University Medical Center, an RSA-based technique was developed and used to measure the micromotion of hip-, knee-, and elbow-prostheses. In Fig. 1, a schematic drawing of the
* Corresponding author. Tel.:#(31) 71 526 3935/1911; fax:#(31) 71 526 6801; e-mail: [email protected]. 0021-9290/98/$19.00 ( 1998 Elsevier Science Ltd. All rights reserved. PII: S0021-9290(98)00025-6
stereo X-ray set-up is presented. The positioning of the tubes causes the X-ray beams to intersect at the position of the implant and leads to a resulting radiograph consisting of two halves. Between the patient and the film cassette, a calibration box is positioned in such a way, that the tantalum markers attached to the upper and lower surfaces of the box are clearly visible in both parts of the radiograph. At the lower surface 26 fiducial markers are attached, including two extra markers to identify the calibration box and to specify whether one
492
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
types of endoprostheses using different types of calibration boxes. The aim of this study was to determine the accuracy and robustness of this newly developed software system and to investigate the influence of film digitization on the overall computation errors of the automated RSA procedure.
2. Materials and methods 2.1. Phantom study
Fig. 1. A schematic drawing of the stereo X-ray setup for the exposure of the left and right parts of the radiograph, needed for the RSA procedure. At the lower surface of the calibration box 26 fiducial markers are attached, including two markers to identify the calibration box and to specify whether one deals with the left or right half of the radiograph. At the upper surface twelve control markers are attached.
deals with the left or right half of the radiograph. At the upper surface of the box, 12 control markers are attached. A baseline stereo-radiograph is taken shortly after the implantation of the prosthesis. This allows the assessment of the initial position of the prosthesis with regard to the bone. At every follow-up examination, the position of the prosthesis relative to the bone is re-assessed. Assuming that the markers in the bone are fixed, and taking the positions of these markers as reference points, one can calculate the spatial translation and rotation of the prosthetic components as a function of time. The total time needed to measure and analyze an RSA radiograph is approximately 1 h. A large share of this time is needed for the determination of the two-dimensional coordinates of markers in both images. Until recently, this was done by manually indicating all the required markers on a high-precision mechanical measurement table (Selvik, 1974; Ryd, 1986; Ka¨rrholm, 1989). To facilitate this time consuming, tedious, and subjective procedure with automated image processing techniques, an analytical software package — implemented on a SUN SparcStation 20 — has been developed in our institute. This package (DIRSA — Digital Interactive RSA) allows the automated detection, identification, and matching of all relevant markers in digitized RSA radiographs and the accurate calculation of the two-dimensional positions of those markers. Furthermore, the calculation of three-dimensional marker coordinates and the calculation of micromotion are integrated in the analytical software package. The developed measurement techniques have been applied to several
To validate the performance of the developed algorithms for the automated detection, identification, and matching of the tantalum markers, the automated RSA procedure was applied to nine RSA radiographs of a constructed phantom. The model consisted of a plexiglass cylinder having a diameter of 50 mm. Eight tantalum markers with a diameter of 1 mm were attached to the cylinder’s surface at well-defined positions. Each phantom radiograph was digitized by a Lumisys 100 laser scanner and by a Vidar VXR-12 CCD scanner. The radiographs were digitized with a pixel size of 174 lm, resulting in a matrix of 2048]2600 pixels (150 DPI) with a gray-scale resolution of 8 bits (Vidar scanner) and 12 bits (Lumisys scanner). After digitization, the images were transferred to the SUN workstation for subsequent processing and analysis. To validate the calibration procedure, i.e. the coordinate transformation from the RSA radiograph to the laboratory coordinate system and the calculation of the positions of both X-ray foci, two error measures were computed: (1) the transformation error, expressed as the mean distance between the known and calculated positions of 13 fiducial markers; and (2) the focus error, expressed as the mean distances from the calculated focus position to eight projection lines, which are defined by the known positions of the control markers and the transformed positions of the projections of those markers. Furthermore, the rigid-body matrices for all phantom radiographs were calculated. A rigid-body matrix contains all distances between cylinder markers in one radiograph. Two occurring distances were selected to validate the automated RSA procedure. The selected distances were 48.00$0.03 and 20.00$0.03 mm, respectively. Both distances occurred four times in each radiograph. The results — calibration errors and rigid body analysis — obtained with both film scanners were compared with the results obtained with the manual procedure. Finally, the relative micromotion of the computed rigid bodies with respect to the position obtained with the measurement table was calculated. Since only rotations about the central axis of the cylinder were applied, the calculated micromotion was expressed in helical axis parameters.
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
2.2. Patient study To assess the in vivo reproducibility of the manual and automated RSA measurements, both procedures were also applied to the stereo-radiographs of 18 patients with an Interax total knee prosthesis. At least three noncollinear tantalum markers were inserted into the tibia (Aronson et al., 1974). The tibial component of the prosthesis was marked with three attached stainless-steel markers with a diameter of 2 mm. The diameter of the tantalum calibration markers was 1 mm, corresponding to approximately 6 pixels in the digitized images. In Fig. 2, a typical stereo-radiograph of the Interax total knee prosthesis is shown. In total, 30 patient radiographs were measured twice with the manual measurement procedure. For the digital procedure, the radiographs were scanned twice with a spatial resolution of 150 DPI and an 8-bit gray scale resolution. Based on our conclusions drawn out of the phantom study, the patient radiographs were digitized with the Vidar VXR-12 only. 2.3. Detection of markers To detect the markers, an extension and improvement of the circle-finding concept described by Duda and Hart (1972) has been developed. Several modifications of the conventional Hough transform to recover circles and ellipses, to increase the computational speed and to decrease memory usage have been described in the literature (Han et al., 1994; Kimme et al., 1975; McKenzie and
493
Protheroe, 1990; Pao et al., 1993; Yip et al., 1992). The application of the Hough transform for the detection of circular objects in a digital image usually requires a three-dimensional array of accumulators. In the Duda and Hart circle-finding algorithm, the array is indexed by three parameters specifying the location and size of a circle. Since the markers are positioned at different object-film distances, the magnification factor is not constant. The differences, however, are only small. Therefore, we could consider the radius of the circles to be approximately known, i.e. lying in the range 1.0— 1.2 mm. The Reason why we reduced the Hough featurespace to a two-dimensional array indexed by the x- and y-coordinate of the position of the circle. The image to which the circle-finding algorithm is applied, is the result of a digital edge enhancement (Sobel operator) applied to the raw digitized data (Duda and Hart, 1973), followed by a pixel thresholding for noise reduction. As a result the input image consists primarily of digitized lines and curves. To determine the positions of the markers with a higher accuracy, the detection of the circular objects as described above is followed by a least-squares error fitting procedure. By fitting a paraboloid to the gray value profile of the markers, the initial estimate of the detected centerpoint is refined to sub-pixel accuracy. The projections of the prosthetic markers, in our case having a diameter of 2 mm, are detected by the same procedure — using a different circle diameter — as used for the detection of the bone and the calibration markers.
Fig. 2. A digitized RSA stereo-radiograph of the Interax total knee prosthesis. A large number of calibration markers, bone markers, and prosthetic markers are clearly visible in the radiograph.
494
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
2.4. Identification of calibration markers From the detected tantalum makers, the calibration markers are extracted and identified. In Fig. 3, the positions of the calibration markers are clearly visible. The first step of the developed identification algorithm involves the detection of the largest rectangle with four markers as corner points. The four corner points as well as all markers lying on the edges or edge extensions of the detected rectangle are labeled as a group, while allowing a certain geometric distortion. In the next step, all individual markers in the detected group must be labeled as upper or lower calibration marker. This is done by matching the known layout of the upper and lower surface of the calibration box to the group of markers. If a model fits, each marker will be labeled as a specific calibration marker. Each pair of markers with the smallest distance is considered to be a correct marker-
model pair. Markers without a model-partner are removed from the group. For a model-marker without a partner, a candidate will be searched among the detected markers that do not belong to the group. In a final step, markers that are incorrectly identified as calibration markers are detected and removed. The described labeling procedure is repeated to label the other plane of the calibration box. The labeling of the detected markers results in three groups: (1) fiducial markers in the lower plane of the calibration box, (2) control markers in the upper plane of the calibration box, and (3) unidentified markers. 2.5. Identification of bone markers and prosthestic markers The unidentified markers mentioned in the previous section are considered to be bone markers. An important step in the RSA procedure is the matching of the bone
Fig. 3. An example of the screen layout of the newly developed software package DIRSA. Detected markers and connections between markers are presented graphically with colored circles and lines. The magnification window, used to correct for the position and label of identified markers, is clearly visible in a popup window.
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
markers and the prosthetic markers in the left and right parts of the radiograph. This matching is carried out by calculating a matrix with elements (i, j) representing the distance between the projection line from the left focus to marker i (in the left radiograph) and the projection line from the right focus to marker j (in the right radiograph). The position halfway the shortest connecting line-piece can be considered as the three-dimensional position of the corresponding marker. The right and left bone marker projection can be matched by selecting the elements that are a row minimum as well as a column minimum. Thereafter, the three-dimensional positions of the markers are calculated. 2.6. Interactive corrections All detection, identification and matching results, obtained with this newly developed RSA software package,
495
can be corrected interactively by the user. All markers can be re-labeled interactively. Several interactive tools are integrated in DIRSA to allow the user to make corrections if the computed results are not satisfactory. Detected and identified markers can be repositioned or re-labeled with the use of a magnification window with a magnification factor up to 64 times (see Fig. 3). Detected markers can be deleted and new markers can be added. The automated matching of the markers in the left and right half of the radiograph can be changed interactively by the user by clicking on the connections in the image (see Fig. 3) and subsequently deleting them or reconnecting them to specific markers. Also, new connections can be added. After an interactive correction, all results are recalculated immediately and automatically. These interactive features are a major advantage of the developed software package.
Fig. 4. An example of the graphical presentation of results computed by DIRSA. All intermediate results are shown in a text window. Furthermore, the calculated rotations and translations in three dimensions are presented in a graphical format.
496
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
2.7. Calculation of micromotion A patient database has been integrated in DIRSA to store relevant patient information and to save the results obtained from the analysis of the RSA radiographs. To calculate the micromotion of prosthetic components, a series of analyzed examinations from one patient can be selected from the database. The software automatically matches the rigid-body matrices in subsequent RSA radiographs. We developed a new algorithm to compare the distances in both matrices and to exchange marker labels to produce two identically labeled matrices. Since markers are not always stably lodged in bone and measurement errors can occur, we allowed a maximum rigid body error of 0.3 mm. After the matching of the rigid bodies of bone and prosthesis, the micromotion of the prosthetic components relative to the bone is calculated. Results can be shown in text format when two examinations are analyzed or in graphical format when a complete series is analyzed. Per analysis two graphs are produced. One graph with the measured translations in the x-, y- and z-direction and one graph with the measured rotations around the x-, y- and z-axis. In Fig. 4, an example of this graphical output is shown.
3. Results 3.1. Phantom study All calibration markers in the phantom radiographs were automatically detected and correctly labeled. All cylinder markers were also detected. In two radiographs, however, the connection between two markers had to be interactively corrected by the observer. The reason was that the computed rigid body matrix contained a series of nearly equal distances, due to the symmetric geometry of the constructed model. In Table 1, the calibration errors of the phantom radiographs are listed for the manual procedure as well as for the automated procedure carried out with the two film scanners. The Lumisys scanner we used showed a severe geometric deformation of the scan raster. The pixel size turned out to be a function of the scanning position. This distortion fluctuated in time and was Table 1 The transformation errors and foci errors obtained from nine phantom radiographs Method
Transformation error (mm)
Focus error (mm)
Manual Lumisys Vidar
0.03 0.51 0.13
0.70 6.66 1.30
therefore difficult to correct for. Due to the resulting inaccuracy in the detection of the calibration box markers, the automated determination of the positions of both X-ray foci showed a considerable error — the crossing-line error was about ten times larger — compared to the manual procedure and compared to the results obtained with the Vidar scanner. In Table 2, the calculated distances between cylinder markers are shown. The results show that the random errors associated with both the automated and manual procedures are very similar. With the automated procedure carried out with the Lumisys scanner, however, a small bias (3%) in distance was found. This is probably caused by the error in the focus position: a systematic difference of about 5% was found between the manually calculated z-coordinate of the X-ray focus position and the automatically computed coordinate. With the Vidar scanner this systematic error was absent and there was no significant difference in the calculated foci positions with regard to the manual procedure. In Table 3, the average motion of the model between subsequent exposures calculated by the automated RSA procedure was compared with the corresponding results of the manual procedure. The motion was expressed in two helical axis parameters, i.e. the rotation about and the translation along the longitudinal axis. Again, the Vidar scanner shows better results than the Lumisys scanner. Based on these experiments, we decided to use the Vidar scanner for the in vivo experiment and also for our clinical studies. 3.2. Patient study Approximately, 85% of all calibration markers and bone markers in the patient radiographs were automatically detected and correctly matched. Table 2 The distances between cylinder markers in the phantom radiographs measured with different techniques
Manual Lumisys Vidar
D (mm) 20
D (mm) 48
20.03$0.02 19.42$0.03 19.96$0.03
48.08$0.02 47.78$0.03 48.03$0.01
Abbreviations: D "distance between markers (actual value 20 20.00 mm); D "distance between markers (actual value 48.00 mm). 48 Table 3 The relative motion of the phantom object in helical axis parameters, calculated from the difference between the manually obtained position and the position obtained by the two automated RSA procedures
Lumisys Vidar
Translation (mm)
Rotation (deg)
0.32$0.17 0.17$0.05
0.18$0.11 0.12$0.11
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498 Table 4 The mean calibration errors and mean rigid-body errors of the manual and automated RSA procedures (n"30) Calibration errors (mm)
Rigid-body errors (mm)
Method
Transform.
Focus
Tibia
Prosthesis
Manual Automated
0.09 0.13
1.19 1.86
0.05 0.05
0.06 0.07
Table 5 The reproducibility of the manual and automated RSA procedures expressed as the standard deviation obtained from repeated measurements in 30 patient radiographs Translation (mm)
Rotation (deg)
Method
D x
D y
D z
/ x
/ y
/ z
Manual Automated
0.03 0.03
0.08 0.04
0.13 0.11
0.26 0.24
0.27 0.23
0.05 0.08
497
ematic analysis — included. With the current version of the software on a SUN SparcStation 20, total processing time is less than 10 min per radiograph (interactive correction procedures included). The software package contains a number of tools to correct for the position and the label of detected markers and to correct the connections between specific markers. Without interaction the total analysis time of one RSA stereo radiograph is even smaller: about 4 min. This processing speed is more or less independent from the total number of markers used. Since 95% of the processing time is spent on the Hough transforms, the detection of e.g. 100 markers is as fast as the detection of just 10 markers. With the increasing speed of computers, the amount of time needed to analyse RSA radiographs with DIRSA will certainly decrease in the near future. At the moment, the software package is being ported to the Windows NT-platform.
4. Discussion and conclusion Abbreviations: D "Translation along the x-axis; D "Translation x y along the y-axis; D "translation along the z-axis; / "Rotation z x about the x-axis; / "Rotation about the y-axis; / "Rotation about y z the x-axis.
The calibration errors and the rigid-body errors, i.e. the difference distances between the markers in the bone or at the prosthesis between two observations, have been summarized in Table 4. The results show that the calibration errors of the automated procedure are comparable with the corresponding errors of the manual procedure. The rigid-body errors also had comparable mean values for both techniques: 0.05 mm for the tibia and 0.06 mm for the prosthesis. The reproducibility of both procedures was assessed by calculating the average standard deviation of the micromotion of the tibial component (which should be zero) between repeatedly measured radiographs. In Table 5, the results are expressed as translations along and rotations about the x-, y- and z-axis. The results show that the reproducibility of the automated procedure is slightly better than the reproducibility of the manual procedure. Due to the orientation of the markers at the prosthesis, i.e. in a plane that is parallel with the radiographic film, the standard deviation of the translation along the z-axis (perpendicular to the film) is the largest. This is also the reason why the deviations of the rotations about the x- and the y-axis are larger than the deviation of the rotation about the z-axis. 3.3. Processing speed The time required for the manual detection and labeling of markers is approximately 1 h for each radiograph, the RSA calculations — rigid-body position and kin-
We developed an analytical software package for the automated detection and identification of markers in RSA radiographs, the computation of the three-dimensional positions of these markers, and for the calculation of the micromotion of endoprostheses. The algorithm for the detection of circular markers has proven to be robust and accurate. False positive detections, occurring in the first step of the detection algorithm, were automatically removed in subsequent steps. Only markers with very low contrast could not be detected automatically, but these were the same markers that caused problems with the manual method. The possibility to add markers ‘manually’ is integrated in our software package. Based on our results, we could say that our detection and identification procedures are very robust. To our opinion, the amount of contrast in the acquired radiograph is a critical factor. When the film contains very dark areas or has a very non-uniform optical density distribution, the software (and also the human eye) can have problems to find the markers, especially the calibration markers at the edge of the radiograph. Usually, one or two markers out of the 13 fiducial markers were not detected. Not all markers of the calibration cage, however, have to be detected to compute accurate migration results. A minimum of seven fiducial markers and four control markers is needed to get good calibration results. The developed detection method is not sensitive to the presence of wires, screws, holes, pins etc. in the radiographs. Only if artificial objects are nearly circular and have the same size as implanted or attached markers, they are detected as candidate markers. After the identification procedure, however, these false positive detections are automatically deleted from the marker list.
498
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
Markers can be implanted very close to each other, without causing serious problems. If two markers overlap more that 50%, however, in most cases only one marker will be detected. In that case, a correction by the user (displace the marker and add a second one) is necessary. The method still works if markers are partly covered by metal or other markers. The local amplitude of the Hough transform depends on the amount of points (gradients) on the edge of a marker. So a marker that is overlapped for 50% by metal will still be detected but half as strong as a complete circle. More than 50% overlap may cause detection problems. The radius of the markers to be detected is a predefined parameter of the Hough transform. In fact, we do not define a constant radius, but a small range (1.0—1.2 mm). So the algorithm uses some inherent ‘averaging’. This causes the algorithm to detect all markers with a radius lying in the predefined range. Using the Vidar CCD scanner, there was no significant difference in calculated focus position compared with the manual RSA procedure. Furthermore, the reproducibility of the automated RSA measurements was comparable to the reproducibility of the manual approach. In conclusion, a new and widely applicable, computerassisted technique has become available to detect, identify and match markers in RSA radiographs and to assess the micromotion of endoprostheses.This new technique will be used in our clinic for our hip, knee, and elbow studies.
References Aronson, A.S., Holst, L., Selvik G., 1974. An instrument for insertion of radiopaque markers. Radiology 113, 733—734. Duda, R.O., Hart, P.E., 1972. Use of the Hough Transformation to detect lines and curves in pictures. Comm. ACM, 11—15. Duda, R.O., Hart, P.E., 1973. Pattern Classifications and Scene Analysis. Wiley, New York. Glasbey, C.A., Horgan, G.W., Hitchcock, D., 1994. A Note on the grey-scale response and sampling properties of a desktop scanner. Pattern Recognition Lett. 15, 705—711. Han, J.H., Ko´czy, L.T., Poston, T., 1994. Fuzzy Hough transform. Pattern Recognition Lett. 15, 649—658. Ka¨rrholm, J., 1989. Roentgen stereophotogrammetry. Review of orthopaedic applications. Acta Orthop. Scand. 60, 491—503. Kimme, C., Ballard, D., Sklansky, J., 1975. Finding circles by an array of accumulators. Commu. ACM 18, 120—122. McKenzie, D.S., Protheroe, S.R., 1990. Curve description using the inverse Hough Transform. Pattern Recognition 23, 283—290. Pao, D., Li, H.F., Jayakumar, R., 1993. A decomposable parameter space for the detection of ellipses. Pattern Recognition Lett. 14, 951—958. Ryd, L., 1986. Micromotion in knee arthroplasty. A Roentgen stereophotogrammetric analysis of tibial component fixation. Acta Orthop. Scand. 57(Suppl.) 220. Selvik, G., 1974. Roentgen Stereophotogrammetry. A method for the study of the kinematics of the skeletal system. PhD. Thesis; reprint Acta Orthop. Scand. 60(Suppl.) 232. So¨derkvist, I., 1990. Some numerical methods for kinematical analysis. PhD. Thesis, University of Umea, UMINF-186.90, ISSN-0348—0542. Yip, R.K.K., Tam, P.K.S., Leung, D.N.S., 1992. Modification of Hough transform for circles and ellipses detection using a 2-dimensional array. Pattern Recognition 25, 1007—1022.
Technical Note
Fast and accurate automated measurements in digitized stereophotogrammetric radiographs Henri A. Vrooman!,*, Edward R. Valstar", Gert-Jan Brand!, Dennis R. Admiraal!, Piet M. Rozing", Johan H.C. Reiber! ! Laboratory for Clinical and Experimental Image Processing, Division of Image processing, Department of Radiology, Leiden University Medical Center, P.O. Box 9600, NL-2300 RC Leiden, The Netherlands " Department of Orthopaedics, Leiden University Medical Center, The Netherlands Received in final form 16 January 1998
Abstract Until recently, Roentgen Stereophotogrammetric Analysis (RSA) required the manual definition of all markers using a highresolution measurement table. To automate this tedious and time-consuming process and to eliminate observer variabilities, an analytical software package has been developed and validated for the detection, identification, and matching of markers in RSA radiographs. The digital analysis procedure consisted of the following steps: (1) the detection of markers using a variant of the Hough circle-finder technique; (2) the identification and labeling of the detected markers; (3) the reconstruction of the three-dimensional position of the bone markers and the prosthetic markers; and (4) the computation of micromotion. To assess the influence of film digitization, the measurements obtained from nine phantom radiographs using two different film scanners were compared with the results obtained by manual processing. All markers in the phantom radiographs were automatically detected and correctly labeled. The best results were obtained with a Vidar VXR-12 CCD scanner, for which the measurement errors were comparable to the errors associated with the manual approach. To assess the in vivo reproducibility, 30 patient radiographs were analyzed twice with the manual as well as with the automated procedure. Approximately, 85% of all calibration markers and bone markers were automatically detected and correctly matched. The calibration errors and the rigid-body errors show that the accuracy of the automated procedure is comparable to the accuracy of the manual procedure. The rigid-body errors had comparable mean values for both techniques: 0.05 mm for the tibia and 0.06 mm for the prosthesis. The reproducibility of the automated procedure showed to be slightly better than that of the manual procedure. The maximum errors in the computed translation and rotation of the tibial component were 0.11 mm and 0.24°, compared to 0.13 mm and 0.27° for the manual RSA procedure. The total processing time is less than 10 min per radiograph, including interactive corrections, compared to approximately 1 h for the manual approach. In conclusion, a new and widely applicable, computer-assisted technique has become available to detect, identify, and match markers in RSA radiographs and to assess the micromotion of endoprostheses. This new technique will be used in our clinic for our hip, knee, and elbow studies. ( 1998 Elsevier Science Ltd. All rights reserved. Keywords: Endoprostheses; Micromotion; Roentgen stereophotogrammetric analysis; Film digitization; Digital image processing
1. Introduction At the Department of Orthopaedics, Leiden University Medical Center, an RSA-based technique was developed and used to measure the micromotion of hip-, knee-, and elbow-prostheses. In Fig. 1, a schematic drawing of the
* Corresponding author. Tel.:#(31) 71 526 3935/1911; fax:#(31) 71 526 6801; e-mail: [email protected]. 0021-9290/98/$19.00 ( 1998 Elsevier Science Ltd. All rights reserved. PII: S0021-9290(98)00025-6
stereo X-ray set-up is presented. The positioning of the tubes causes the X-ray beams to intersect at the position of the implant and leads to a resulting radiograph consisting of two halves. Between the patient and the film cassette, a calibration box is positioned in such a way, that the tantalum markers attached to the upper and lower surfaces of the box are clearly visible in both parts of the radiograph. At the lower surface 26 fiducial markers are attached, including two extra markers to identify the calibration box and to specify whether one
492
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
types of endoprostheses using different types of calibration boxes. The aim of this study was to determine the accuracy and robustness of this newly developed software system and to investigate the influence of film digitization on the overall computation errors of the automated RSA procedure.
2. Materials and methods 2.1. Phantom study
Fig. 1. A schematic drawing of the stereo X-ray setup for the exposure of the left and right parts of the radiograph, needed for the RSA procedure. At the lower surface of the calibration box 26 fiducial markers are attached, including two markers to identify the calibration box and to specify whether one deals with the left or right half of the radiograph. At the upper surface twelve control markers are attached.
deals with the left or right half of the radiograph. At the upper surface of the box, 12 control markers are attached. A baseline stereo-radiograph is taken shortly after the implantation of the prosthesis. This allows the assessment of the initial position of the prosthesis with regard to the bone. At every follow-up examination, the position of the prosthesis relative to the bone is re-assessed. Assuming that the markers in the bone are fixed, and taking the positions of these markers as reference points, one can calculate the spatial translation and rotation of the prosthetic components as a function of time. The total time needed to measure and analyze an RSA radiograph is approximately 1 h. A large share of this time is needed for the determination of the two-dimensional coordinates of markers in both images. Until recently, this was done by manually indicating all the required markers on a high-precision mechanical measurement table (Selvik, 1974; Ryd, 1986; Ka¨rrholm, 1989). To facilitate this time consuming, tedious, and subjective procedure with automated image processing techniques, an analytical software package — implemented on a SUN SparcStation 20 — has been developed in our institute. This package (DIRSA — Digital Interactive RSA) allows the automated detection, identification, and matching of all relevant markers in digitized RSA radiographs and the accurate calculation of the two-dimensional positions of those markers. Furthermore, the calculation of three-dimensional marker coordinates and the calculation of micromotion are integrated in the analytical software package. The developed measurement techniques have been applied to several
To validate the performance of the developed algorithms for the automated detection, identification, and matching of the tantalum markers, the automated RSA procedure was applied to nine RSA radiographs of a constructed phantom. The model consisted of a plexiglass cylinder having a diameter of 50 mm. Eight tantalum markers with a diameter of 1 mm were attached to the cylinder’s surface at well-defined positions. Each phantom radiograph was digitized by a Lumisys 100 laser scanner and by a Vidar VXR-12 CCD scanner. The radiographs were digitized with a pixel size of 174 lm, resulting in a matrix of 2048]2600 pixels (150 DPI) with a gray-scale resolution of 8 bits (Vidar scanner) and 12 bits (Lumisys scanner). After digitization, the images were transferred to the SUN workstation for subsequent processing and analysis. To validate the calibration procedure, i.e. the coordinate transformation from the RSA radiograph to the laboratory coordinate system and the calculation of the positions of both X-ray foci, two error measures were computed: (1) the transformation error, expressed as the mean distance between the known and calculated positions of 13 fiducial markers; and (2) the focus error, expressed as the mean distances from the calculated focus position to eight projection lines, which are defined by the known positions of the control markers and the transformed positions of the projections of those markers. Furthermore, the rigid-body matrices for all phantom radiographs were calculated. A rigid-body matrix contains all distances between cylinder markers in one radiograph. Two occurring distances were selected to validate the automated RSA procedure. The selected distances were 48.00$0.03 and 20.00$0.03 mm, respectively. Both distances occurred four times in each radiograph. The results — calibration errors and rigid body analysis — obtained with both film scanners were compared with the results obtained with the manual procedure. Finally, the relative micromotion of the computed rigid bodies with respect to the position obtained with the measurement table was calculated. Since only rotations about the central axis of the cylinder were applied, the calculated micromotion was expressed in helical axis parameters.
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
2.2. Patient study To assess the in vivo reproducibility of the manual and automated RSA measurements, both procedures were also applied to the stereo-radiographs of 18 patients with an Interax total knee prosthesis. At least three noncollinear tantalum markers were inserted into the tibia (Aronson et al., 1974). The tibial component of the prosthesis was marked with three attached stainless-steel markers with a diameter of 2 mm. The diameter of the tantalum calibration markers was 1 mm, corresponding to approximately 6 pixels in the digitized images. In Fig. 2, a typical stereo-radiograph of the Interax total knee prosthesis is shown. In total, 30 patient radiographs were measured twice with the manual measurement procedure. For the digital procedure, the radiographs were scanned twice with a spatial resolution of 150 DPI and an 8-bit gray scale resolution. Based on our conclusions drawn out of the phantom study, the patient radiographs were digitized with the Vidar VXR-12 only. 2.3. Detection of markers To detect the markers, an extension and improvement of the circle-finding concept described by Duda and Hart (1972) has been developed. Several modifications of the conventional Hough transform to recover circles and ellipses, to increase the computational speed and to decrease memory usage have been described in the literature (Han et al., 1994; Kimme et al., 1975; McKenzie and
493
Protheroe, 1990; Pao et al., 1993; Yip et al., 1992). The application of the Hough transform for the detection of circular objects in a digital image usually requires a three-dimensional array of accumulators. In the Duda and Hart circle-finding algorithm, the array is indexed by three parameters specifying the location and size of a circle. Since the markers are positioned at different object-film distances, the magnification factor is not constant. The differences, however, are only small. Therefore, we could consider the radius of the circles to be approximately known, i.e. lying in the range 1.0— 1.2 mm. The Reason why we reduced the Hough featurespace to a two-dimensional array indexed by the x- and y-coordinate of the position of the circle. The image to which the circle-finding algorithm is applied, is the result of a digital edge enhancement (Sobel operator) applied to the raw digitized data (Duda and Hart, 1973), followed by a pixel thresholding for noise reduction. As a result the input image consists primarily of digitized lines and curves. To determine the positions of the markers with a higher accuracy, the detection of the circular objects as described above is followed by a least-squares error fitting procedure. By fitting a paraboloid to the gray value profile of the markers, the initial estimate of the detected centerpoint is refined to sub-pixel accuracy. The projections of the prosthetic markers, in our case having a diameter of 2 mm, are detected by the same procedure — using a different circle diameter — as used for the detection of the bone and the calibration markers.
Fig. 2. A digitized RSA stereo-radiograph of the Interax total knee prosthesis. A large number of calibration markers, bone markers, and prosthetic markers are clearly visible in the radiograph.
494
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
2.4. Identification of calibration markers From the detected tantalum makers, the calibration markers are extracted and identified. In Fig. 3, the positions of the calibration markers are clearly visible. The first step of the developed identification algorithm involves the detection of the largest rectangle with four markers as corner points. The four corner points as well as all markers lying on the edges or edge extensions of the detected rectangle are labeled as a group, while allowing a certain geometric distortion. In the next step, all individual markers in the detected group must be labeled as upper or lower calibration marker. This is done by matching the known layout of the upper and lower surface of the calibration box to the group of markers. If a model fits, each marker will be labeled as a specific calibration marker. Each pair of markers with the smallest distance is considered to be a correct marker-
model pair. Markers without a model-partner are removed from the group. For a model-marker without a partner, a candidate will be searched among the detected markers that do not belong to the group. In a final step, markers that are incorrectly identified as calibration markers are detected and removed. The described labeling procedure is repeated to label the other plane of the calibration box. The labeling of the detected markers results in three groups: (1) fiducial markers in the lower plane of the calibration box, (2) control markers in the upper plane of the calibration box, and (3) unidentified markers. 2.5. Identification of bone markers and prosthestic markers The unidentified markers mentioned in the previous section are considered to be bone markers. An important step in the RSA procedure is the matching of the bone
Fig. 3. An example of the screen layout of the newly developed software package DIRSA. Detected markers and connections between markers are presented graphically with colored circles and lines. The magnification window, used to correct for the position and label of identified markers, is clearly visible in a popup window.
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
markers and the prosthetic markers in the left and right parts of the radiograph. This matching is carried out by calculating a matrix with elements (i, j) representing the distance between the projection line from the left focus to marker i (in the left radiograph) and the projection line from the right focus to marker j (in the right radiograph). The position halfway the shortest connecting line-piece can be considered as the three-dimensional position of the corresponding marker. The right and left bone marker projection can be matched by selecting the elements that are a row minimum as well as a column minimum. Thereafter, the three-dimensional positions of the markers are calculated. 2.6. Interactive corrections All detection, identification and matching results, obtained with this newly developed RSA software package,
495
can be corrected interactively by the user. All markers can be re-labeled interactively. Several interactive tools are integrated in DIRSA to allow the user to make corrections if the computed results are not satisfactory. Detected and identified markers can be repositioned or re-labeled with the use of a magnification window with a magnification factor up to 64 times (see Fig. 3). Detected markers can be deleted and new markers can be added. The automated matching of the markers in the left and right half of the radiograph can be changed interactively by the user by clicking on the connections in the image (see Fig. 3) and subsequently deleting them or reconnecting them to specific markers. Also, new connections can be added. After an interactive correction, all results are recalculated immediately and automatically. These interactive features are a major advantage of the developed software package.
Fig. 4. An example of the graphical presentation of results computed by DIRSA. All intermediate results are shown in a text window. Furthermore, the calculated rotations and translations in three dimensions are presented in a graphical format.
496
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
2.7. Calculation of micromotion A patient database has been integrated in DIRSA to store relevant patient information and to save the results obtained from the analysis of the RSA radiographs. To calculate the micromotion of prosthetic components, a series of analyzed examinations from one patient can be selected from the database. The software automatically matches the rigid-body matrices in subsequent RSA radiographs. We developed a new algorithm to compare the distances in both matrices and to exchange marker labels to produce two identically labeled matrices. Since markers are not always stably lodged in bone and measurement errors can occur, we allowed a maximum rigid body error of 0.3 mm. After the matching of the rigid bodies of bone and prosthesis, the micromotion of the prosthetic components relative to the bone is calculated. Results can be shown in text format when two examinations are analyzed or in graphical format when a complete series is analyzed. Per analysis two graphs are produced. One graph with the measured translations in the x-, y- and z-direction and one graph with the measured rotations around the x-, y- and z-axis. In Fig. 4, an example of this graphical output is shown.
3. Results 3.1. Phantom study All calibration markers in the phantom radiographs were automatically detected and correctly labeled. All cylinder markers were also detected. In two radiographs, however, the connection between two markers had to be interactively corrected by the observer. The reason was that the computed rigid body matrix contained a series of nearly equal distances, due to the symmetric geometry of the constructed model. In Table 1, the calibration errors of the phantom radiographs are listed for the manual procedure as well as for the automated procedure carried out with the two film scanners. The Lumisys scanner we used showed a severe geometric deformation of the scan raster. The pixel size turned out to be a function of the scanning position. This distortion fluctuated in time and was Table 1 The transformation errors and foci errors obtained from nine phantom radiographs Method
Transformation error (mm)
Focus error (mm)
Manual Lumisys Vidar
0.03 0.51 0.13
0.70 6.66 1.30
therefore difficult to correct for. Due to the resulting inaccuracy in the detection of the calibration box markers, the automated determination of the positions of both X-ray foci showed a considerable error — the crossing-line error was about ten times larger — compared to the manual procedure and compared to the results obtained with the Vidar scanner. In Table 2, the calculated distances between cylinder markers are shown. The results show that the random errors associated with both the automated and manual procedures are very similar. With the automated procedure carried out with the Lumisys scanner, however, a small bias (3%) in distance was found. This is probably caused by the error in the focus position: a systematic difference of about 5% was found between the manually calculated z-coordinate of the X-ray focus position and the automatically computed coordinate. With the Vidar scanner this systematic error was absent and there was no significant difference in the calculated foci positions with regard to the manual procedure. In Table 3, the average motion of the model between subsequent exposures calculated by the automated RSA procedure was compared with the corresponding results of the manual procedure. The motion was expressed in two helical axis parameters, i.e. the rotation about and the translation along the longitudinal axis. Again, the Vidar scanner shows better results than the Lumisys scanner. Based on these experiments, we decided to use the Vidar scanner for the in vivo experiment and also for our clinical studies. 3.2. Patient study Approximately, 85% of all calibration markers and bone markers in the patient radiographs were automatically detected and correctly matched. Table 2 The distances between cylinder markers in the phantom radiographs measured with different techniques
Manual Lumisys Vidar
D (mm) 20
D (mm) 48
20.03$0.02 19.42$0.03 19.96$0.03
48.08$0.02 47.78$0.03 48.03$0.01
Abbreviations: D "distance between markers (actual value 20 20.00 mm); D "distance between markers (actual value 48.00 mm). 48 Table 3 The relative motion of the phantom object in helical axis parameters, calculated from the difference between the manually obtained position and the position obtained by the two automated RSA procedures
Lumisys Vidar
Translation (mm)
Rotation (deg)
0.32$0.17 0.17$0.05
0.18$0.11 0.12$0.11
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498 Table 4 The mean calibration errors and mean rigid-body errors of the manual and automated RSA procedures (n"30) Calibration errors (mm)
Rigid-body errors (mm)
Method
Transform.
Focus
Tibia
Prosthesis
Manual Automated
0.09 0.13
1.19 1.86
0.05 0.05
0.06 0.07
Table 5 The reproducibility of the manual and automated RSA procedures expressed as the standard deviation obtained from repeated measurements in 30 patient radiographs Translation (mm)
Rotation (deg)
Method
D x
D y
D z
/ x
/ y
/ z
Manual Automated
0.03 0.03
0.08 0.04
0.13 0.11
0.26 0.24
0.27 0.23
0.05 0.08
497
ematic analysis — included. With the current version of the software on a SUN SparcStation 20, total processing time is less than 10 min per radiograph (interactive correction procedures included). The software package contains a number of tools to correct for the position and the label of detected markers and to correct the connections between specific markers. Without interaction the total analysis time of one RSA stereo radiograph is even smaller: about 4 min. This processing speed is more or less independent from the total number of markers used. Since 95% of the processing time is spent on the Hough transforms, the detection of e.g. 100 markers is as fast as the detection of just 10 markers. With the increasing speed of computers, the amount of time needed to analyse RSA radiographs with DIRSA will certainly decrease in the near future. At the moment, the software package is being ported to the Windows NT-platform.
4. Discussion and conclusion Abbreviations: D "Translation along the x-axis; D "Translation x y along the y-axis; D "translation along the z-axis; / "Rotation z x about the x-axis; / "Rotation about the y-axis; / "Rotation about y z the x-axis.
The calibration errors and the rigid-body errors, i.e. the difference distances between the markers in the bone or at the prosthesis between two observations, have been summarized in Table 4. The results show that the calibration errors of the automated procedure are comparable with the corresponding errors of the manual procedure. The rigid-body errors also had comparable mean values for both techniques: 0.05 mm for the tibia and 0.06 mm for the prosthesis. The reproducibility of both procedures was assessed by calculating the average standard deviation of the micromotion of the tibial component (which should be zero) between repeatedly measured radiographs. In Table 5, the results are expressed as translations along and rotations about the x-, y- and z-axis. The results show that the reproducibility of the automated procedure is slightly better than the reproducibility of the manual procedure. Due to the orientation of the markers at the prosthesis, i.e. in a plane that is parallel with the radiographic film, the standard deviation of the translation along the z-axis (perpendicular to the film) is the largest. This is also the reason why the deviations of the rotations about the x- and the y-axis are larger than the deviation of the rotation about the z-axis. 3.3. Processing speed The time required for the manual detection and labeling of markers is approximately 1 h for each radiograph, the RSA calculations — rigid-body position and kin-
We developed an analytical software package for the automated detection and identification of markers in RSA radiographs, the computation of the three-dimensional positions of these markers, and for the calculation of the micromotion of endoprostheses. The algorithm for the detection of circular markers has proven to be robust and accurate. False positive detections, occurring in the first step of the detection algorithm, were automatically removed in subsequent steps. Only markers with very low contrast could not be detected automatically, but these were the same markers that caused problems with the manual method. The possibility to add markers ‘manually’ is integrated in our software package. Based on our results, we could say that our detection and identification procedures are very robust. To our opinion, the amount of contrast in the acquired radiograph is a critical factor. When the film contains very dark areas or has a very non-uniform optical density distribution, the software (and also the human eye) can have problems to find the markers, especially the calibration markers at the edge of the radiograph. Usually, one or two markers out of the 13 fiducial markers were not detected. Not all markers of the calibration cage, however, have to be detected to compute accurate migration results. A minimum of seven fiducial markers and four control markers is needed to get good calibration results. The developed detection method is not sensitive to the presence of wires, screws, holes, pins etc. in the radiographs. Only if artificial objects are nearly circular and have the same size as implanted or attached markers, they are detected as candidate markers. After the identification procedure, however, these false positive detections are automatically deleted from the marker list.
498
H.A. Vrooman et al. / Journal of Biomechanics 31 (1998) 491—498
Markers can be implanted very close to each other, without causing serious problems. If two markers overlap more that 50%, however, in most cases only one marker will be detected. In that case, a correction by the user (displace the marker and add a second one) is necessary. The method still works if markers are partly covered by metal or other markers. The local amplitude of the Hough transform depends on the amount of points (gradients) on the edge of a marker. So a marker that is overlapped for 50% by metal will still be detected but half as strong as a complete circle. More than 50% overlap may cause detection problems. The radius of the markers to be detected is a predefined parameter of the Hough transform. In fact, we do not define a constant radius, but a small range (1.0—1.2 mm). So the algorithm uses some inherent ‘averaging’. This causes the algorithm to detect all markers with a radius lying in the predefined range. Using the Vidar CCD scanner, there was no significant difference in calculated focus position compared with the manual RSA procedure. Furthermore, the reproducibility of the automated RSA measurements was comparable to the reproducibility of the manual approach. In conclusion, a new and widely applicable, computerassisted technique has become available to detect, identify and match markers in RSA radiographs and to assess the micromotion of endoprostheses.This new technique will be used in our clinic for our hip, knee, and elbow studies.
References Aronson, A.S., Holst, L., Selvik G., 1974. An instrument for insertion of radiopaque markers. Radiology 113, 733—734. Duda, R.O., Hart, P.E., 1972. Use of the Hough Transformation to detect lines and curves in pictures. Comm. ACM, 11—15. Duda, R.O., Hart, P.E., 1973. Pattern Classifications and Scene Analysis. Wiley, New York. Glasbey, C.A., Horgan, G.W., Hitchcock, D., 1994. A Note on the grey-scale response and sampling properties of a desktop scanner. Pattern Recognition Lett. 15, 705—711. Han, J.H., Ko´czy, L.T., Poston, T., 1994. Fuzzy Hough transform. Pattern Recognition Lett. 15, 649—658. Ka¨rrholm, J., 1989. Roentgen stereophotogrammetry. Review of orthopaedic applications. Acta Orthop. Scand. 60, 491—503. Kimme, C., Ballard, D., Sklansky, J., 1975. Finding circles by an array of accumulators. Commu. ACM 18, 120—122. McKenzie, D.S., Protheroe, S.R., 1990. Curve description using the inverse Hough Transform. Pattern Recognition 23, 283—290. Pao, D., Li, H.F., Jayakumar, R., 1993. A decomposable parameter space for the detection of ellipses. Pattern Recognition Lett. 14, 951—958. Ryd, L., 1986. Micromotion in knee arthroplasty. A Roentgen stereophotogrammetric analysis of tibial component fixation. Acta Orthop. Scand. 57(Suppl.) 220. Selvik, G., 1974. Roentgen Stereophotogrammetry. A method for the study of the kinematics of the skeletal system. PhD. Thesis; reprint Acta Orthop. Scand. 60(Suppl.) 232. So¨derkvist, I., 1990. Some numerical methods for kinematical analysis. PhD. Thesis, University of Umea, UMINF-186.90, ISSN-0348—0542. Yip, R.K.K., Tam, P.K.S., Leung, D.N.S., 1992. Modification of Hough transform for circles and ellipses detection using a 2-dimensional array. Pattern Recognition 25, 1007—1022.