Intensity-Based Nonoverlapping Area Registration Supporting “Drop-Outs” in Terms of Model-Based Radiostereometric Analysis

A model-based radiostereometric analysis (MBRSA) is a method for precise measurement of prosthesis migration, which does not require marking the implant with tantalum beads. Instead, the prosthesis pose is typically recovered using a feature-based 2D-3D registration of its virtual model into a stereo pair of radiographs. In this study, we evaluate a novel intensity-based formulation of previously published nonoverlapping area (NOA) approach. The registration is capable of performing with both binary radiographic segmentations and nonsegmented X-ray images. In contrast with the feature-based version, it is capable of dealing with unreliable parts of prosthesis. As the straightforward formulation allows efficient acceleration using modern graphics adapters, it is possible to involve precise high-poly virtual models. Moreover, in case of binary segmentations, the nonoverlapping area is simply interpretable and useful for indicating the accuracy of the registration outcome. In silico and phantom evaluations were performed using a cementless Zweymüller femoral stem and its reverse engineered (RE) model. For initial pose estimates with difference from the ground-truth limited to ±4 mm and ±4°, respectively, the mean absolute translational error was not higher than 0.042 ± 0.035 mm. The error in rotation around the proximodistal axis was 0.181 ± 0.265°, and the error for the remaining axes was not higher than 0.035 ± 0.037°.


Introduction
Radiostereometric analysis (RSA), introduced by Selvik [1,2], is an established method for an accurate measurement of prosthesis mechanical stability, indicated in particular in cases of total joint arthroplasty. e analysis is used for measuring micromotion between the prosthesis and the surrounding bone. Due to its high precision, it allows to reveal a potential failure of the implant fixation at early stages, when the prosthesis migration is not recognizable in plain radiographs, nor clinical symptoms occur [3]. e conventional radiostereometric analysis depends on two sets of tantalum beads. e first set of markers is attached to the prosthesis, while the second set of beads is injected directly into a bone surrounding the implant. e position of each marker in three-dimensional space is obtained using a triangulation from a stereo pair of radiographs. Commonly, a patient undergoes several following-up radiographic examinations during a certain time period after the arthroplasty [4]. A potential failure of the prosthesis fixation is observed when the relative pose between the two sets of markers differs between the individual examinations.
However, the attachment of tantalum beads to the implant raises several potential issues. In radiographs, the prosthesis may occlude the attached beads, the marked implants are more expensive, and the strength of the prosthesis may be negatively affected. To overcome these difficulties, model-based radiostereometric analysis (MBRSA) has been proposed by Valstar et al. [5]. Instead of attaching the beads, the implant pose is recovered by 2D-3D registration of its virtual model into a stereo pair of radiographs. Several studies have revealed that the model-based radiostereometry reaches lower but acceptable accuracy in comparison with the conventional approach [6][7][8].
Registration methods used in radiostereometry are typically feature-based, exploiting edges detected in radiographs and a prosthesis outline obtained from the virtual model. Valstar et al. [5] proposed an approach based on nonoverlapping area (NOA) minimization, which required a complete outline of the prosthesis to be obtained from the radiographs. e major drawback of the method was an inability to handle unreliable parts of the detected outline, as there were significantly large dimensional differences between the actual prosthesis and its CAD model involved in the phantom evaluation. A following-up study, proposed by Kaptein et al. [9], enhanced the accuracy by using reversed engineering (RE) models of prosthesis instead of CAD models provided by a manufacturer and by improving the registration method to handle the unreliable parts of detected contours, referred to as "drop-outs," which may be caused by metallic objects that are not included in the virtual model, such as bone screws. e registration was based on minimization of contour difference, which can be in contrast with the original nonoverlapping area evaluated locally, and the unreliable parts of the contour, selected by the user, may be simply omitted from the registration. e minimization of contour difference was in broader principle adopted by many subsequent studies [6,7,10,11].
In this study, we propose an intensity-based radiostereometric method, reviving the idea of nonoverlapping area. In contrast with [5], the proposed registration allows to evaluate the nonoverlapping area locally. Consequently, the contribution of this revisited method is the ability to handle the drop-outs and unreliable parts of the prosthesis captured in radiographs. As the contour detection and a feature matching are not required by the intensity-based registration, the computation is much more simple in comparison with the previously published approaches. erefore, the method is straightforward for efficient acceleration using graphics adapters. e study presents in silico and phantom evaluations of the proposed approach.  Figure 2.

Materials
Final positions of the tantalum markers with respect to the prosthesis were calculated, as their locations inside the phantom were defined by the CAD model, used for manufacturing the box. e Carestream Directview DR 9500 System was exploited for sequential capturing of digital radiographs (DR). e phantom assembly was inserted into a biplanar calibration cage filled with 36 tantalum beads. We used direct linear transformation (DLT) [12] for the radiograph calibration, as proposed by Choo and Oxland [13], instead of the traditional fiducial and control planes (FCP) approach [2]. e assembly was rotated approximately 45°around the prosthesis proximodistal axis to prevent occlusions of the phantom markers by the implant. e complete experimental setup is shown in Figure 3.
For the phantom study, 8 radiographs were captured from each anterior-posterior and lateral views. e pose of the calibration box within the imaging system was varied among the individual acquirements. e radiographs were enhanced using an intensity curve adjustment and histogram equalization. Upon the radiographs, a set of 64 stereo pairs was constructed, and an example stereo pair of radiographs is shown in Figure 4. Randomly chosen 32 pairs were exploited for a precise refinement of the mutual pose between the phantom and the implant, and the remaining half was used for the evaluation.

Methods
e proposed registration is suitable for usage with both binary segmentations and enhanced nonsegmented radiographs.
3.1. Binary Images. As the metallic implants are highly radiopaque, the segmentation is performed by thresholding the enhanced radiographs; pixels representing the prosthesis are set to 1. A coarse initial estimate of the prosthesis pose must be provided by the user. During the registration, binary digitally reconstructed radiographs (DRR) are rendered from the prosthesis model. Following [5], the nonoverlapping area is defined as the area that the segmentations of real and calculated radiographs do not have in common. e size of the area is equal to the count of different pixels between the real and virtual segmentations. Since the segmentations contain only binary values, the count is computed by summing squares of the pixel differences (PD): where p � (R, T) is a vector formed by a rotation and translation of the prosthesis model in the space of stereo radiographs. To eliminate different radiograph resolutions or perspective scaling, it is convenient to express the nonoverlapping area size in a relative form as NOA(p)/ (NOA(p) + C(p)), where C(p) is a count of overlapping pixels. e metric is schematically depicted in Figure 5. e minimization of the nonoverlapping area for anterior-posterior (AP) and lateral (LAT) views is formulated as nonlinear leastsquares (NLS) problem:

Nonsegmented Radiographs.
Due to the significant radiopacity, it is assumed that the metallic prostheses are objects with the highest contrast in radiographs, exceeding the brightness of the surrounding bone, soft tissues, or eventual cement layer, which makes the segmentation rather a straightforward task. On the other hand, a precise segmentation may demand some additional user interaction, and consequently, to decrease the amount of required user efforts, it is convenient to perform the registration using directly the nonsegmented radiographs. In case of the proposed intensity-based registration, radiographs are preprocessed using the histogram equalization. After the preprocessing, pixels representing the prosthesis reach approximately the maximum value of the image intensity range. e digitally reconstructed radiographs contain only two intensity levels. Following the radiopacity assumption, the virtual model is rendered with the highest contrast, while the background pixels are set to the lowest intensity. As the intensity-based nonoverlapping area registration is formulated as a least-squares problem, it is clearly suitable for usage with gray-scale images. With respect to the high prosthesis radiopacity and the least-squares formulation of the registration, the equalized radiographs may be considered as    a probabilistic approximation of the prosthesis segmentation. However, in this case, the sum of squared differences does not correspond to the exact size of the nonoverlapping area, in contrast with the registration involving only the binary segmentations.

Handling Drop-Outs.
Drop-outs are especially related with metallic objects that are not a part of the prosthesis virtual model, but which are present in radiographs and occlude certain parts of the implant. In case of hip prosthesis, the ball head attached to the femoral stem may be occluded by a metallic acetabular implant. In this case, a user must roughly select the area, where a boundary of the prosthesis, corresponding to the virtual model, is unclear. e situation is schematically illustrated in Figure 6. Consequently, the drop-out areas are discarded from both input X-ray images and digitally reconstructed radiographs; hence, they do not affect the registration accuracy. e drop-outs are supported by both segmentation-based and intensitybased registrations. However, they were not supported by the original contour-based approach [5], as it required a complete and precise outline of the prosthesis to be extracted from the input radiographs.

Optimization Scheme.
During the registration, 6 degrees of freedom (DOF) of the model pose are optimized using Levenberg-Marquardt numerical solver [14]. As the Levenberg-Marquardt optimization is gradient based, an evaluation of Jacobian matrix J F is required during each iteration [15]. e matrix contains partial derivatives of pixel differences with respect to the pose parameters: e number of rows in J F matrix is given by the total count of pixels in both anterior-posterior and lateral images, and the number of columns is equal to the count of degrees of freedom. As it is not possible to evaluate the Jacobian matrix using a closed-form solution, we use a central difference approximation: where p ∈ p is a certain pose parameter and ε is a difference spacing. To increase both capture range and accuracy at the same time, the registration is divided into five subsequent optimizations where the coarse-to-fine strategy is applied on the different spacing ε. Stages with ε equal to 1e 1 , 1e 0 , 1e −1 , 1e −2 , and 1e −3 millimeters or degrees, respectively, were used in the study. To speed up the registration and lower the memory requirements, only regions of interest were cropped from radiographs to form the pixel differences vector, based on bounding boxes of the implant segmentations.
To prevent an undesirable cropping of the nonoverlapping areas, the bounding boxes were enlarged by certain margins. Due to pixel-wise formulation of the registration, places containing drop-outs, selected by the user, were simply discarded from the registration.

In Silico Evaluation.
As eventual segmentation errors may negatively affect the 2D-3D registration [16], the aim of the in silico evaluation was to investigate the intensity-based approach accuracy itself, without external influences. ree data sets containing one hundred virtual stereo radiographs of the implant, differing in resolution, were created with pixel spacings set to 0.5, 0.35, and 0.143 mm. e initial poses were generated randomly with uniform distribution, and the maximal translational and rotational errors were limited to ± 5 mm and ± 5°, respectively. Table 1 shows mean values and standard deviations of absolute pose errors together with corresponding nonoverlapping area size, number of iterations, and processing time. A relation between accuracy and pixel spacing is shown using box plots in Figure 7. e accuracy obviously increases with the radiograph resolution, as the registration is able to perform more iterations. However, the rising count of iterations together with increasing length of the pixel differences (PD) vector yields into a trade-off between the registration accuracy and the processing time.

Phantom Evaluation.
e model-based radiostereometric analysis monitors possible changes in relative pose between the bone and the joint replacement among certain time periods. e prosthesis pose is recovered by the registration of its virtual model into a stereo pair of radiographs, while the pose of the bone is obtained using a tantalum beads placed inside the bone. e tantalum beads are inserted using commercially available injectors, provided for instance by Tilly Medical Products AB or RSA Biomedical Suppliers. A threedimensional pose of the bone markers is easily obtained from   the stereo radiographs by triangulation. In consequence, the prosthesis migration is measured with respect to the set of markers injected into the bone. During the phantom study, for the accuracy evaluation purposes, the ground-truth pose of the implant within the space of stereo radiographs was determined using ten tantalum markers inside the phantom Plexiglas box, as the relative pose between the phantom and the attached prosthesis was known. e registration was evaluated for both binary and nonsegmented radiographic images, with and without user selected drop-outs. A sample stereo pair containing drop-outs, chosen from the evaluation data set, is shown in Figure 8. e results of accuracy evaluations with and without drop-outs are shown in Tables 1 and 2, respectively, and their comparison is visualized in Figure 9. To investigate a relation between the capture range and the registration accuracy, the evaluations were performed for different limitations of maximal errors in initial poses, revealing a gentle decrease of accuracy and higher number of iterations with rising initial pose error.
e results show that the estimation of proximodistal rotation reaches the lowest accuracy in comparison with other pose parameters. e accuracy of the rotation around the  y-axis would be increased by involving a third radiographic image taken in proximodistal projection, allowing the registration to minimize a nonoverlapping area even in the xz plane. However, in a real clinical environment, it is not possible to capture a radiographic image from such projection.
Generally, the recovery of the prosthesis pose using its virtual model is possible due to sufficient asymmetry of the implant, leading to unique projections of the model [5,9]. erefore, dropping the ball head out from the radiographs, a significantly asymmetric part of the prosthesis, which may be on the other hand in real situation occluded by a metallic acetabular implant, has rather slight but still recognizable influence on the registration accuracy. As the ball head is the most proximal and the most medial part of the model, a slight decrease of the accuracy can be seen mainly in the translation along the proximodistal and z-axis. e method also performs for binary segmentations with slightly higher accuracy than for enhanced radiographs. e registration pipeline was implemented using Qt Toolkit 5.8.0 and compiled with MSVC 2013 64-bit. To speedup the registration, the rendering part of the pipeline was accelerated using OpenGL 4.3. e evaluations were performed using a Microsoft Windows 8.1 64-bit desktop machine equipped with an Intel Core i5-6500 CPU processor, an NVidia 980 GTX Ti 6GB graphics adapter, and a 24 GB DDR4 SDRAM memory module.

Discussion and Conclusion
In contrast with 2D-3D registration methods exploiting contour difference minimization, the nonoverlapping area does not require feature matching between detected and virtual contours, which is a computationally demanding and error-prone task. Considering theoretical aspects, determining matches between the contours is in principle an ill-posed problem. Strictly speaking, there are no actual correspondences between the detected and calculated contours until the ground-truth pose of the model is recovered. In other poses, the virtual contour captures different places of the prosthesis than the edges detected in radiographs. We therefore suggest that the nonoverlapping area has a stronger theoretical basis than the contour difference registration.
A computation of the intensity-based nonoverlapping area is more straightforward in comparison with the original feature-based formulation. In the feature-based case, the area was evaluated using nontrivial procedure based on horizontal directions of both detected and virtual contours [5]. On the contrary, the intensity-based variant exploits plain pixel differences between radiographs and virtual segmentations obtained from the prosthesis model.
As the OpenGL acceleration was focused only on the part of the pipeline, data transfers between a graphical and operational memory were a cause of a performance bottleneck. ere is an opportunity for further significant acceleration by implementing the rest of the registration pipeline using the OpenGL compute shader programs, eliminating the memory transfers and exploiting parallelization of the Levenberg-Marquardt algorithm. We believe the shift from the feature to intensity-based variant is possible due to rapid progress of hardware performance, as the intensity-based registration feasibility depends on usage of modern hardware resources.
Due to efficient graphics hardware and intensity-based formulation, it is possible to involve complete high-poly RE models without decimating the mesh, in contrast with studies presented by Kaptein et al. and Seehaus et al. [9,17]. e registration accuracy is comparable with previously published feature-based approaches, according to the summary presented by Syu et al. [11]. However, the comparison is rather tentative, as the accuracy depends on the shape of involved implants [6] and on the type of imaging system. An important contribution of the intensity-based revision is the ability to handle the dropouts, which are useful for dealing with components that are not a fixed part of the prosthesis model. We also suggest that the relative size of the nonoverlapping area is a simply interpretable metric, useful for indicating the resulting accuracy of the registration.

Conflicts of Interest
e authors declare no conflicts of interest.