Model-based Roentgen stereophotogrammetry of orthopaedic implants

Journal of Biomechanics 34 (2001) 715–722 Model-based Roentgen stereophotogrammetry of orthopaedic implants E.R. Valstara,*, F.W. de Jonga, H.A. Vroomanb, P.M. Rozinga, J.H.C. Reiberb b a Department of Orthopaedics, Leiden University Medical Center, P.O. Box 9600, 2300 RC Leiden, The Netherlands Division of Image Processing, Department of Radiology, Leiden University Medical Center, P.O. Box 9600, 2300 RC Leiden, The Netherlands Accepted 14 February 2001 Abstract Attaching tantalum markers to prostheses for Roentgen stereophotogrammetry (RSA) may be difficult and is sometimes even impossible. In this study, a model-based RSA method that avoids the attachment of markers to prostheses is presented and validated. This model-based RSA method uses a triangulated surface model of the implant. A projected contour of this model is calculated and this calculated model contour is matched onto the detected contour of the actual implant in the RSA radiograph. The difference between the two contours is minimized by variation of the position and orientation of the model. When a minimal difference between the contours is found, an optimal position and orientation of the model has been obtained. The method was validated by means of a phantom experiment. Three prosthesis components were used in this experiment: the femoral and tibial component of an Interax total knee prosthesis (Stryker Howmedica Osteonics Corp., Rutherfort, USA) and the femoral component of a Profix total knee prosthesis (Smith & Nephew, Memphis, USA). For the prosthesis components used in this study, the accuracy of the model-based method is lower than the accuracy of traditional RSA. For the Interax femoral and tibial components, significant dimensional tolerances were found that were probably caused by the casting process and manual polishing of the components surfaces. The largest standard deviation for any translation was 0.19 mm and for any rotation it was 0.528. For the Profix femoral component that had no large dimensional tolerances, the largest standard deviation for any translation was 0.22 mm and for any rotation it was 0.228. From this study we may conclude that the accuracy of the current model-based RSA method is sensitive to dimensional tolerances of the implant. Research is now being conducted to make model-based RSA less sensitive to dimensional tolerances and thereby improving its accuracy. # 2001 Elsevier Science Ltd. All rights reserved. Keywords: Roentgen stereophotogrammetry; Three-dimensional pose estimation; Surface models; Micromotion 1. Introduction Roentgen stereophotogrammetric analysis (RSA) is an accurate measurement technique to assess micromotion of implants with respect to the surrounding bone (Selvik, 1989). In RSA, the three-dimensional position and orientation of objects is determined by the reconstruction of the three-dimensional position of well-defined markers. For this purpose, tantalum markers are used that are inserted into the bone and are either attached to or inserted into the implant. However, marking of implants may be difficult and is sometimes even impossible. Furthermore, marking of implants is an expensive procedure and in some *Corresponding author. Tel.:+31-71-5262975; 5266743. E-mail address: [email protected] (E.R. Valstar). fax:+31-71- countries it is only allowed by the regulatory bodies after extensive testing and comprehensive documentation. For some implants, like metal-backed cups in total hip arthroplasty and femoral components in total knee arthroplasty, the metal of the implant often obscures the attached markers when RSA radiographs are taken. RSA studies of these implants with attached markers are only possible when care is taken in positioning the patient during radiography. Because of these difficulties, only one clinical RSA study of femoral components in total knee arthroplasty has been performed (Nilsson et al., 1995). In contrast, many clinical RSA studies of tibial components in total knee arthroplasty have been conducted (Overview in K.arrholm, 1989). Several attempts have been made to perform RSA studies without attaching markers. For several clinical RSA studies of hip stems, the head of the prosthesis has . been used as a marker (Onsten et al., 1995; K.arrholm, 0021-9290/01/$ - see front matter # 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 1 - 9 2 9 0 ( 0 1 ) 0 0 0 2 8 - 8 716 E.R. Valstar et al. / Journal of Biomechanics 34 (2001) 715–722 1989, K.arrholm et al., 1997). For polyethylene cups with a metal ring, the ring has been used to obtain the cup’s position (Snorrason and K.arrholm, 1990). Furthermore, for hemispherical metal-backed cups the position as well as the orientation could be assessed by using the projection of the hemispherical part and the projection of the base circle (Valstar et al., 1997). Others have used the shape of hip stems to obtain the position and the orientation of the stem (Turner-Smith and Bulstrode, 1993; Valstar, 1996, 2001). All of these techniques used basic geometrical shapes}circles, spheres, and straight lines or single well defined landmarks to define the position and orientation of the implant. These techniques cannot be used for total knee prostheses and other more complex shaped implants since the shape of these implants cannot be described by these basic geometrical shapes. A model-based RSA method has been developed to overcome the above mentioned problems. With this method the three-dimensional position and orientation of complex shaped prosthesis components is assessed without the use of markers. The method is based on matching of the detected contour of an implant (that is inserted in a patient) with the calculated projected contour of a three-dimensional model of the same implant (De Jong, 1997). Similar techniques have been used for other applications: the determination of the position of vertebrae (Laval!ee and Szeliski, 1997) and the assessment of the position and orientation of total knee prostheses from single focus fluoroscopic images (Banks and Hodge, 1996; Walker et al., 1996; Zuffi et al., 1999). In this study, the model-based RSA method is presented and the accuracy of the method is tested by an in vitro experiment with a phantom. 2. Material and method 2.1. Material A phantom study was carried out with a femoral component and a tibial component of an Interax total knee prosthesis (Stryker Howmedica Osteonics Corp., Rutherfort, USA) and a femoral component of a Profix total knee prosthesis (Smith & Nephew, Memphis, USA). The phantom was a Plexiglas cylinder with a diameter of 40 mm with 12 1 mm spherical tantalum markers embedded in its surface. These 12 tantalum markers were used to define a local coordinate system. For each experiment, one of the prosthesis components was rigidly attached to the base plane of this cylinder. The phantom was placed in an RSA set-up that consisted of two synchronized Roentgen tubes that were positioned at approximately 1.5 m above a film cassette. Each Roentgen tube was directed at one half of the film under an angle with the vertical of approximately 208. A Plexiglas calibration box with 1 mm tantalum markers was positioned underneath the Roentgen table. This calibration box defined the three-dimensional laboratory coordinate system and was used to accurately calculate the foci positions. The x-axis of this coordinate system was directed in the medio-lateral direction, the y-axis was directed in the caudo-cranial direction, and the z-axis was directed in the posteroanterior direction. Within this RSA set-up, the phantom was positioned in seven successive poses and at each pose an RSA radiograph was taken. Subsequent radiographs were analyzed and the relative position and orientation of the implant with respect to the cylinder was assessed. By using its attached tantalum markers, the position and orientation of the cylinder could be assessed very accurately (Valstar et al., 1997, 2000). Since the prosthesis and the cylinder were rigidly connected, their actual relative motion was zero. When comparing the relative position and orientation of the prosthesis and the cylinder between two successive stereoradiographs within a series, changes in position and orientation of the implant relative to the cylinder could be assessed. Since the actual relative motion was zero, these changes indicated the error of the model-based RSA method. 2.2. Method Radiography is based on the central projection of objects on radiographic film by X-rays that originate from an X-ray focus. A central projection is described mathematically by projection parameters. In RSA, these projection parameters are assessed by a calibration procedure. The first step in the model-based RSA method is the detection of the contour of the implant in the radiograph. Thereafter, the projected contour of a threedimensional prosthesis model is calculated using the projection parameters that were calculated during the calibration procedure. Finally, the pose of the implant is calculated by minimizing the difference between the detected contour and the calculated model contour. Since the projection of an object with sufficient asymmetry is unique, the position and orientation of the implant have been assessed. However, implants have dimensional tolerances and will therefore deviate from the three-dimensional model. Together with lacking flatness of the X-ray film and measurement errors, these deviations will result in a detected contour and a calculated model contour that will not fit exactly. Nevertheless, a minimal difference between the two contours corresponds with an optimal approximation of the position and orientation of the implant. E.R. Valstar et al. / Journal of Biomechanics 34 (2001) 715–722 2.2.1. Model creation There are many ways to describe the three-dimensional shape of an object. For this application we have chosen to describe the object as a triangulated surface model. The surface of these models is represented by a mesh that is composed of a large number of triangles (elements) and nodes. A node is a geometrical location defined by its coordinates. The sequence of the nodes in an element determines the direction of the normal vector of that element. To assure the connectivity of the elements, the normal vectors of all elements should have the same direction with respect to the object: either inward or outward. The triangulated surface models were created with MSC/PATRAN (MacNeal–Schwendler GmbH, . Munchen, Germany). CAD (Computer Aided Design)models of the implants, which were provided by the manufacturers of the implants, was used as input for MSC/PATRAN. In Fig. 1 such a CAD-model and the resulting triangulated surface model are shown. The accuracy of the triangulated surface model depends on the accuracy of the input that was used to create the model, on the number of triangles used, and on the distribution of the triangles. The higher the curvature of the object, the smaller the local triangle-size should be. There is no absolute criterion for the required number of elements of the model, but in this study approximately 5000 elements were used. Details in the prosthesis that never will be part of the outer contour of the model, such as blind holes or threading are redundant. To increase the speed of the algorithm, these details were not modeled. 2.2.2. Contour detection The contours of the implant in the radiograph are detected by means of the Canny operator (Canny, 1986). After thresholding of the gradient strength, a set of binary morphologic operations, such as dilation and skeletonization, is applied to extract the major closed contours. Except for the outer contour a projected Fig. 1. The solid model of the Interax femoral component together with its meshed representation. 717 implant can also have several inner contours, such as holes in the implant. The direction of the outer contour is defined clockwise, inner contours are defined counterclockwise. 2.2.3. Model contour creation The elements of the triangulated surface model can be represented by a support plane: n
x ¼ d; jnj ¼ 1 in which n is the normal vector of the element, x is an arbitrary point, and d is a scalar. A given element is defined visible when the focus is situated in front of the element’s support plane and invisible when the focus is situated in or behind the element’s support plane. The following equations hold: ðn
xf  dÞ > 0 ðvisibleÞ; ðn
xf  dÞ40 ðinvisibleÞ; in which xf is the focus position. As the elements describe a closed surface, each edge has exactly two elements attached to it. The edges in the model that are connected to a visible and an invisible element are termed contour-edges. These contour-edges always form one or more closed chains in three-dimensional space. The result of the projection of such a three-dimensional contour is a set of closed chains of two-dimensional nodes. In contrast to the detected contour, the calculated model contour consists of closed chains that may intersect (Fig. 2). 2.2.4. Non overlapping area The non overlapping area (NOA) is defined as the area that the detected contour and the calculated model contour do not have in common. To determine the NOA, both contours are divided into intervals indicated by the x-coordinates of the nodes in the calculated model contour (Fig. 3). The number of intervals thus equals the number of contour nodes minus 1. On each Fig. 2. A detected image contour and a calculated model contour in a detail of an RSA radiograph of the Interax tibial component. 718 E.R. Valstar et al. / Journal of Biomechanics 34 (2001) 715–722 function are summed: NOAtotal ¼ NOAleft þNOAright : Fig. 3. Determination of area code of the NOA. In order to simplify the explanation in this figure the inner contours of the calculated model contour have been omitted. The arrows indicate the direction of the contours, the vertical lines indicate an area between two fictive nodes, and the numbers indicate the area code. The shaded area, with area 1, represents the NOA. interval the direction of the contour is known. Parts of the contour with a positive x-direction are assigned direction value 1, parts with a negative x-direction are assigned direction value 1. An area code that has the same value as the direction value is assigned to the area above each contour part. Areas beneath contour parts receive area code zero. The area codes have to be calculated for both the detected contour and the calculated model contour. By summing the area codes of both areas a resulting area code is obtained. If the summed area code is –1 or 1, the area is added to the NOA. 2.2.5. Minimization of the non overlapping area In order to calculate the three-dimensional position and orientation of the implant, the NOA has to be minimized. The minimization of the NOA is carried out by an optimization procedure In this procedure an objective function, F, is used to find the position and orientation of the model resulting in a minimal NOA. In RSA two projections of an object are used: one projection in the left image half and one projection in the right image half. For both image halves an object function may be defined as: NOAleft ¼ FðModelpose ; Model; Focusleft ; Contourleft Þ; NOAright ¼ FðModelpose ; Model; Focusright ; Contourright Þ: In these two objective functions, the model pose is the only unknown, since the model of the prosthesis is given, the focus position has been assessed during the calibration procedure, and the contours have been detected in the image before the optimization procedure started. In order to obtain a better stability and accuracy of the optimization procedure, the two objective The model pose consists of six parameters}three position parameters and three orientation parameters }and is highly non-linear. In each step of the minimization procedure, the current position and orientation of the model is used to calculate a new model contour. The NOA is derived by using the procedure as described in the previous paragraph. For minimization of the non-linear objective function a minimization scheme called Feasible Sequential Quadratic Programming was used (FSQP; Lawrence and Tits, 1996; Panier and Tits, 1993). This scheme was implemented by using the C-library CSFQP, which may be downloaded from: http://www.isr.umd.edu/Labs/CACSE/FSQP/ fsqp.html The model-based RSA method was incorporated in the DIRSA software that has been described in Vrooman et al. (1998). This software was designed to run on a SUN-workstation. In this experiment a SUN SparcStation 20 (Sun Microsystems, Inc., Palo Alto, Ca, USA) was used. 2.3. Results In Fig. 4, one of the radiographs from the phantom study of the Interax tibial component is shown. The manually positioned rectangles indicate the region of interest for the contour detection algorithm. The line adjacent to the contour of the tibial component indicates the automatically detected contour of the component. In Fig. 5, a detail of the left image half of Fig. 4 is used to illustrate the optimization procedure. After the actual contour has been detected (Fig. 5a) the optimization procedure starts with an initial estimate of the prosthesis position and orientation (Fig. 5b). The optimization procedure minimizes the NOA and in Fig. 5c an intermediate result is shown. In Fig. 5d the final result of the optimization is presented: a minimal NOA that results in an optimal estimation of the three-dimensional position and orientation of the model. The results of the phantom experiment are summarized in Table 1. The largest mean error for translation was found for the Interax tibial component and was 0.157 mm. Also for rotation the largest mean error was found for the tibial component: 0.3418. Furthermore, the S.D. of the rotations were rather large. The largest S.D. was found for the rotation about the y-axis of the Interax femoral component: 0.5248. For the tibial component the S.D. was the largest for the rotation about the same coordinate axis: 0.4278. For the Profix femoral component, the mean errors of the translation parameters were smaller than the mean E.R. Valstar et al. / Journal of Biomechanics 34 (2001) 715–722 719 Fig. 4. A stereoradiograph (projected on a single film) of the cylinder with attached Interax tibial component. errors found for both Interax components. For the rotation parameters the mean errors as well as the S.D. were smaller for the Profix femoral component. The largest standard deviation, 0.2178 was found for the rotation about the z-axis. The rather large standard deviations in the parameters of the Interax femoral component were probably caused by large dimensional tolerances on the medial side of the component resulting in a mismatch between the projected contour and the calculated model contour (Fig. 6a). An explanation for this mismatch is that, during manufacturing, implants are being polished by hand after they have been cast in a mold. The casting process and the manual polishing process do not provide constant dimension tolerances along the prosthesis’ entire surface. In areas that require a high accuracy, like the patella gliding surface of the femoral component, the dimension tolerances will be much smaller than in regions that do not require such a high accuracy, like the edges of the implant. However, in this case the size of the mismatch was so large, about 1.5 mm, that it could indicate that an inaccurate CADmodel has been used. In the Interax tibial component, poor dimensional tolerances also resulted in rather large errors. The tibial component consists of two parts that are fitted into each other during surgery. The conical shaped hole in the tibial baseplate that fitted the cruciform part had poor dimensional tolerances. As a result, the direction of this hole deviated approximately 18 from the ideal perpendicular orientation with respect to the tibial baseplate (Fig. 6b). Although for the Profix femoral component no large dimensional tolerances were found, some smaller mismatches were found in this component (Fig. 6c). The edges of the actual prosthesis are rounded, where the edges of the model contour have a more pronounced square shape. The differences between actual implant and model are, however, much smaller than found for both Interax components, and as a result the errors were also smaller. The processing time of a standard RSA radiograph is between 5 and 8 min when tantalum markers are used (Vrooman et al., 1998). This time includes detection of all markers}calibration box markers, bone markers, and prosthesis markers and interactive corrections of intermediate results. The optimization procedure described in this article takes an additional 45 s up to 5 min. The time needed for the optimization depends on the closeness of the first estimate to the final solution. In the worst case, the maximum total analysis time per RSA radiograph was still within 15 min. 3. Discussion and conclusion For two out of three components tested in this study, rather large errors in position and orientation 720 E.R. Valstar et al. / Journal of Biomechanics 34 (2001) 715–722 Table 1 The micromotion of the Interax femoral and tibial components and the Profix femoral component relative to the cylinder markers. Translations are labeled x, y, and z (in mm) and rotations are labeled Rx , Ry , and Rz (in8; n=7) x y z Rx Ry Rz Interax femur Mean 0.096 0.045 0.098 0.214 0.258 0.112 S.D. 0.190 0.080 0.165 0.353 0.524 0.146 Interax tibia Mean 0.157 0.012 0.022 0.341 0.064 0.250 S.D. 0.075 0.070 0.182 0.213 0.427 0.195 Profix femur Mean 0.047 0.016 0.023 0.020 0.075 0.074 S.D. 0.203 0.133 0.221 0.145 0.173 0.217 Fig. 5. The optimization process with the Interax tibial component illustrated by the projected contour in the left half of the RSA image. (a) After the region of interest has been determined by the observer, the contour of the prosthesis is automatically detected by means of the Canny operator. (b) The first estimation of the position and orientation of the prosthesis model is projected in the radiograph. (c) An intermediate result of the optimization procedure: the overlap of both contours is increasing. (d) The final results: of the optimization procedure: an optimum overlap of both contours has been obtained. were found that were probably caused by large dimensional tolerances in certain areas of the component’s surfaces. The results of these two components were not as good as reported for RSA studies that used prostheses with attached markers (Ryd, 1986; Nilsson et al., 1991, 1995). The standard deviations of repeated measurements that were reported in these RSA studies ranged between 0.032 and 0.15 mm for translations and 0.07 and 0.138 for rotations, compared to 0.08–0.22 mm for translations and 0.15– 0.528 for rotations found in this study. Especially, for the rotations of the two Interax components large S.D. were observed. However, for one of the three components tested, the Profix femoral component, smaller dimensional differences between the actual prosthesis and the model were observed and the micromotion results were more accurate. The mean values of the micromotion parameters, especially for the rotations, were closer to zero than observed for the Interax components. Also the S.D. for the rotations were smaller than for both Interax components, however, they were still larger than the standard deviations that were observed for prostheses that had attached markers. In order to be an alternative for RSA studies of implants with attached markers, this model-based RSA method should have a better accuracy than the accuracy that now has been reached. This higher accuracy could be obtained by making the method insensitive to dimensional tolerances. This could be done by omitting parts of the detected contour that belong to areas of the implant that have large dimensional tolerances. Another way to increase the accuracy of this model-based RSA method could be the adaptation of the models for each individual implant. In order to obtain the additional data that is necessary for this approach, several points on the surface of the implant will have to be measured by an accurate three-dimensional measuring device. With this data, the triangulated surface model of the implant may be scaled and an individual model for each implant may be created. E.R. Valstar et al. / Journal of Biomechanics 34 (2001) 715–722 721 Fig. 6. Due to large dimensional tolerances a rather large mismatch occurs at the medial side of the Interax femoral component (a) and at the cruciform part of the Interax tibial component (b). The mismatch for the Profix femoral component is much smaller (c) and exists only in the edges of the component. From this study we may conclude that the accuracy of the current model-based RSA method is sensitive to dimensional tolerances of the implant. We aim at improvement of this model-based RSA method so that it will be less sensitive to large dimensional tolerances and that it will provide an accuracy that is comparable to the accuracy of traditional RSA. Currently, these improvements to the model-based RSA method are subject of research that is carried out at our department. References Banks, S.A., Hodge, W.A., 1996. Accurate measurement of threedimensional knee replacement kinematics using single-plane fluoroscopy. IEEE Transactions on Biomedical Engineering 43, 638–649. Canny, A., 1986. Computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence 8, 679– 698. De Jong, F.W., 1997. Mechanical analysis of a hydroxyl-apatite coated knee prosthesis. Thesis, Technical University Eindhoven, The Netherlands, ISBN 90-3860149-2. K.arrholm, J., 1989. Roentgen stereophotogrammetry. Review of orthopaedic applications. Acta Orthopaedica Scandinavica 60 (4), 491–503. K.arrholm, J., Herberts, P., Hultmark, P., Malchau, H., Nivbrant, B., Thanner, J., 1997. Radiostereometry of hip prostheses. Review of methodology and clinical results. Clinical Orthopaedics and Related Research 344, 94–110. Laval!ee, S., Szeliski, R., 1997. Recovering the position and orientation of free-form objects from image contour using 3D distance maps. IEEE Transactions on Pattern Analysis and Machine Intelligence 17, 378–390. 722 E.R. Valstar et al. / Journal of Biomechanics 34 (2001) 715–722 Lawrence, C.T., Tits, A.L., 1996. Nonlinear Equality Constraints in Feasible Sequential Quadratic Programming. Optimization Methods and Software 6, 265–282. Nilsson, K.G., K.arrholm, E.L., Magnusson, P., 1991. Evaluation of micromotion in cemented vs uncemented knee arthroplasty in osteoarthrosis and rheumatoid arthritis: Randomized study using roentgen stereophotogrammetric analysis. Journal of arthroplasty 6 (3), 265–278. Nilsson, K.G., K.arrholm, J., Linder, L., 1995. Femoral component migration in total knee arthroplasty: Randomized study comparing cemented and uncemented fixation of the Miller-Galante design. Journal of Orthopaedic Research 13 (3), 347–356. ( . Onsten, I., Akeson, K., Obrant, K.J., 1995. Migration of the Charnley stem in rheumatoid arthritis and osteoarthritis. A roentgen stereophotogrammetric and histomorphometric study. Journal of Bone and Joint Surgery, British Volume 77 (1), 677–680. Panier, E., Tits, A.L., 1993. On Combining Feasibility, Descent and Superlinear Convergence In Inequality Constrained Optimization. Mathematical Programming 59, 261–276. Ryd, L., 1986. Micromotion in knee arthroplasty: A roentgen stereophotogrammetric analysis of tibial component fixation. Acta Orthopaedica Scandinavica 57 (Suppl 220), 1–80. Selvik, G., 1989. Roentgen stereophotogrammetry: a method for the study of kinematics of the skeletal system. Acta orthopaedica Scandinavica 60 (Suppl 232), 1–51. Snorrason, F., K.arrholm, J., 1990. Early loosening of revision hip arthroplasty. A roentgen stereophotogrammetric analysis. Journal of Arthroplasty 5 (3), 647–651. Turner-Smith, A.R., Bulstrode, C.J., 1993. Stereoradiogrammetry for the prediction of hip replacement survival. In:Turner-Smith, A.R. (Ed.), Micromovement in Orthopaedics, Oxford Medical Engineering Series; vol 10. ISBN 0-19-262306-0:113-123. Valstar, E.R., 1996. Roentgen stereophotogrammetric analysis of the Mallory Head prosthesis without markers. Acta Orthopaedica Scandinavica 67 (Suppl 272), 3. Valstar, E.R., Spoor, C.W., Nelissen, R.G.H.H., Rozing, P.M., 1997. Roentgen stereophotogrammetric analysis of metal-backed cups without attached markers. Journal of Orthopaedic Research 15 (6), 869–873. Valstar, E.R., 2001. Digital Roentgen stereophotogrammetry: Development, validation, and clinical application. Ph.D. Thesis, Leiden University, ISBN 90-9014397-1. Valstar, E.R., Vrooman, H.A., Toksvig-Larsen, S., Ryd, L., Nelissen, R.G.H.H., 2000. Digital automated RSA compared to manually operated RSA. Journal of Biomechanics 33 (12), 1593–1599. Vrooman, H.A., Valstar, E.R., Brand, G.J., Admiraal, D.R., Rozing, P.M., Reiber, J.H.C., 1998. Fast and accurate automated measurements in digitized stereophotogrammetric radiographs. Journal of Biomechanics 31 (5), 491–498. Walker, S.A., Hoff, W., Komistek, R., Dennis, D., 1996. In vivo pose estimation of artificial knee implants using computer vision. Biomedical Sciences Instrumentation 32, 143–150. Zuffi, S., Leardini, A., Catani, F., Fantozzi, S., Cappello, A., 1999. A Model-Based Method for the Reconstruction of Total Knee Replacement Kinematics. IEEE Transactions on Medical Imaging 18, 981–991.