A method to reconstruct patient-specific proximal femur surface models from planar pre-operative radiographs

Medical Engineering & Physics 32 (2010) 1180–1188

Contents lists available at ScienceDirect

Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy

A method to reconstruct patient-specific proximal femur surface models from planar pre-operative radiographs P.E. Galibarov, P.J. Prendergast, A.B. Lennon ∗ Trinity Centre for Bioengineering, School of Engineering, Trinity College Dublin, Dublin 2, Ireland

a r t i c l e

i n f o

Article history: Received 28 May 2010 Received in revised form 20 July 2010 Accepted 17 August 2010 Keywords: Patient-specific 3D reconstruction Image processing Finite element analysis

a b s t r a c t Three-dimensional reconstruction from volumetric medical images (e.g. CT, MRI) is a well-established technology used in patient-specific modelling. However, there are many cases where only 2D (planar) images may be available, e.g. if radiation dose must be limited or if retrospective data is being used from periods when 3D data was not available. This study aims to address such cases by proposing an automated method to create 3D surface models from planar radiographs. The method consists of (i) contour extraction from the radiograph using an Active Contour (Snake) algorithm, (ii) selection of a closest matching 3D model from a library of generic models, and (iii) warping the selected generic model to improve correlation with the extracted contour. This method proved to be fully automated, rapid and robust on a given set of radiographs. Measured mean surface distance error values were low when comparing models reconstructed from matching pairs of CT scans and planar X-rays (2.57–3.74 mm) and within ranges of similar studies. Benefits of the method are that it requires a single radiographic image to perform the surface reconstruction task and it is fully automated. Mechanical simulations of loaded bone with different levels of reconstruction accuracy showed that an error in predicted strain fields grows proportionally to the error level in geometric precision. In conclusion, models generated by the proposed technique are deemed acceptable to perform realistic patient-specific simulations when 3D data sources are unavailable. © 2010 IPEM. Published by Elsevier Ltd. All rights reserved.

1. Introduction In clinical practice medical images are widely used for diagnostic, pre-operative planning and other purposes. Common imaging techniques such as CT, MRI, ultrasound, and conventional radiography are capable of capturing important data for applications, where patient-specific morphology is required. Volumetric images (CT, ␮CT, MRI, 3D ultrasound) are highly desired by researchers who are interested in precise anatomical morphology. However, these images are not always available for research and clinical use. Often just planar radiographs are used for pre-operative planning due to their cost efficiency and lower radiation dose [1]. Due to the lack of 3D data it is more difficult to obtain a precise geometry from such a planar image. Conventional 3D reconstruction from a volumetric image begins with a stack of 2D images (slices) with known spatial dimensions. This process consists of several steps: (i) image

∗ Corresponding author at: Trinity Centre for Bioengineering, School of Engineering, Parsons Building, Trinity College Dublin, Dublin 2, Ireland. Tel.: +353 01 896 2396; fax: +353 01 679 5554. E-mail addresses: [email protected], [email protected] (A.B. Lennon).

processing—enhancing initial images to remove noise and scanning artefacts, (ii) image segmentation—segmenting 2D/3D images into several regions representing objects or, alternatively, a 3D image can be converted directly into a 3D model using voxels (e.g. volume ray casting [2] or finite element models for stress analysis [3]); and, finally, (iii) surface extraction—defining the objects surface representation. An example of surface extraction is the generation of isosurfaces based on voxel intensities, e.g. the Marching Cubes algorithm [4] (closest voxels with a specified Hounsfield value [5] are connected into volumetric polygons to organize a 3D surface.) Alternatively, 2D contours extracted from image slices can be swept into a single surface as has been done by many researchers [6–9]. Many open-source and commercial software packages are available for solving these tasks in an automated and semi-automated manner, e.g. Slicer (Harvard Medical School, Boston, MA, USA), 3D Doctor (Able Software, Lexington, MA, USA), Amira (Visage Imaging Inc., Richmond, Australia), etc. In cases when high quality volumetric data is not available statistical shape modelling (SSM) can be employed to reconstruct a body part. Barratt et al. [10] proposed a technique to register a mean statistical femur model to 3D ultrasound image data. Comparison of the results to corresponding CT-based models showed that they were favourable compared to other existing techniques (average

1350-4533/$ – see front matter © 2010 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2010.08.009

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188

error of 3.1 mm for SSM compared to a range of 2.94–6.88 mm for the other methods [11,12]). Heimann et al. [13] reviewed the application of SSM to different organs. For instance, Fleute et al. [14] used sparse CT images (reduced number of slices) combined with SSM to reconstruct patient-specific femur geometries [14,15]. A study by Shim et al. [15] showed that root mean square errors for reconstructed femur geometry varied between 0.3 and 2.53 mm (subject to the number of slices). An alternative approach to 3D reconstruction using medical imaging data is to use conventional 2D images which do not explicitly contain 3D information, e.g. X-rays, to generate an inverse transform from the projection that created the planar image. As only a single or small number of images may be available, extra knowledge about the shape to be reconstructed is required to create a 3D surface. Most methods within this approach differ in their assumptions and/or algorithm applied to bridging the gap between 2D projection to 3D volume. Bras et al. [16] proposed an approach based on low-dose bi-planar radiographs and non-stereo corresponding contours (NSCC). Zanetti et al. [17] used both single and bi-planar radiographs to morph template femurs to a patientspecific model. An alternative method based on SSM was proposed by Zheng and co workers in a series of studies [18–20]. Non-rigid registration was used to find a set of matched 2D point pairs between features extracted from X-rays and features extracted from a 3D point distribution model. These 2D point pairs were then used to create a set of 3D point pairs using a ray projection algorithm between the 3D model’s silhouette and the extracted image contour. More established 3D registration algorithms could then be employed to achieve the final patient-specific 3D models. Langton et al. [21,22] based their method on a projection of a 3D CT image to create simulated single plane DXA images of bone mineral density and used General Procrustes Analysis of selected landmarks to find and align the 2D image and 3D template model. Thin plate splines were applied to warp the template model. Bredbenner and Nicolella [23] have also used projection of 3D CT data to generate simulated BMD images but used Principal Component Analysis (PCA) to determine a template statistical shape model. Although some of the above studies have used single radiographs to generate their 2D–3D transform, they were either assessed on isolated cadaveric femora, did not consider the femur specifically, or required at least some manual pre-processing of images to separate the femur from other features on the radiograph. Therefore, the assessment of these methods does not account for the difficulty of acquiring and automatically processing femur X-rays in vivo, e.g. overlapping of pelvis and proximal femur on X-ray, and, thus, have to be further tested in realistic pre-operative conditions. Another important issue in using models based on these reconstruction techniques is that no analysis has been performed to date on the effect of geometric inaccuracy on the mechanical parameters predicted by finite element analyses using such models. Such data would be very beneficial to analysts in deciding if a 2D–3D reconstruction technique can be sufficient for their needs in a given study or whether it is necessary to acquire new 3D data to meet their needs. Given the current almost exclusive use of planar radiographs in clinical practice for many procedures, an automated method for generation of femur geometry from a single planar radiograph would undoubtedly be of use in a large number of patient-specific modelling tasks. It is also essential for retrospective studies used to corroborate surgical planning tools [24]. The current study aims to develop a technique to generate extracortical geometry of a patient’s femur from a pre-operative 2D X-ray in an automated manner. Since the lack of 3D data on a planar X-ray will most likely introduce errors to the final geometry an estimate for the error compared to CT-based reconstruction is given. Finally, a conclusion needs to be made as to the suitability of these geometries for patient-specific modelling. In particular, the influence of geomet-

1181

Fig. 1. Method scheme for 3D femur reconstruction from planar radiograph.

ric precision on mechanical responses, such as strain prediction in a finite element analysis, is evaluated. 2. Methods The following strategy was adopted to solve the task of generation of extracortical femur geometry from planar radiographs (see Fig. 1): (i) extract femur contours from a planar radiograph, (ii) select the closest model from a library of generic models, and (iii) warp this model to fit the patient’s femur. In the following the main details are presented in order to maintain a concise presentation of the methods used. A supplemental document containing a more detailed description of the reconstruction methodology accompanies the main article and can be downloaded from the publisher’s website or requested directly from the corresponding author. 2.1. Femur contour extraction algorithm Digital pelvic anterior–posterior (AP) X-rays with known magnification factor are used to reconstruct patient geometries. Often only the proximal femur is visible on pre-operative X-rays. Thus, the present technique targets reconstruction of the proximal part of an individual femur. The algorithm assumes a left femur as input. If reconstruction of a right femur is required the X-ray is reflected and the reconstructed 3D model is later reflected to return to the right hand side. To extract contours of the femur on an X-ray image the following steps are taken: (i) the radiograph is enhanced, (ii) anatomical landmarks are found to construct an initial contour, and (iii) the contour is transformed using an Active Shape algorithm [25,26] to fit the femur. Radiograph enhancement is performed using the Image Processing toolbox, Matlab (The MathWorks, Inc., Natick, MA, USA). First, adaptive contrast enhancement algorithms are applied to adjust initial contrast and brightness. An adaptive low-pass noise-removal filter based on Wiener’s method [27] eliminates noise pixels. Edgedetection methods (Sobel’s method [27] combined with vertical and horizontal gradients) are then applied and the resulting binary image is processed using morphological operations—isolated pixels are removed and small objects are connected. Pixel analysis of the binary image combined with knowledge of femur anatomy [28,29], size, and image orientation enables identification of parts of the femur boundaries. First, a ray-tracing algorithm, in which rays are directed horizontally and record edge crossings, is used to find the distal borders of the femur. According to the table of distributions by Noble et al. [28], minimum extra-

1182

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188

Fig. 2. (a) Schematic femur (adopted from Noble et al. [28]), and landmarks for predicting femoral head centre as in Sugano et al. [29] (red) and (b) landmarks (red) used to reconstruct initial femur contour (purple). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

cortical width of a femur is 20.5 mm and maximum canal width at the lesser trochanter area +20 mm is 63 mm (see Fig. 2a). This information is used to discriminate unlikely candidate points. Iterative neighbourhood pixel analysis reveals lateral and medial borders. Next, all image pixels outside the borders are eliminated as they do not contain required points of the contour. The lateral femur border is used to locate the greater trochanter as the most laterally distant point from the least squares line fitted into the lateral border points. A cropped binary image is used to find contours of the femoral head via a Hough transform [30]. It was suggested by Sugano et al. [29] that a best prediction of the head centre should be based on the midpoint of the femoral neck isthmus. When it is not available the position of the head centre can be estimated from the height of the femoral neck saddle and the distance from the medullary axis to the medial cortex at the 30% level (Fig. 2a). Using these landmarks, along with the known extracortical borders, a convex hull shape can be defined containing the femur. The remaining segments of the extracortical contour are extracted using an Active Contour algorithm [25,26]. This algorithm is an image segmentation technique based on iterative morphing of a 2D closed contour using pixel intensities as external forces. A gradient vector flow (GVF) implementation [25] was chosen due to robustness in case of concave regions. This implementation uses an intensity-based gradient vector field as potential forces, which minimize the energy functional given by Eq. (1).

Eas =

1  1 2



(˛| (s)|2 + ˇ| (s)|2 ) + Eext (s) ds

(1)

0

where (s) ={x(s), y(s)}, s ∈ [0, 1] is the closed contour, and x(s), y(s) are corresponding point coordinates on the curve; ˛ and ˇ are tension and rigidity parameters, respectively and Eext is the energy of external forces. Tension and rigidity parameters were found empir-

ically and equal to 1 and 0.3 respectively. A starting active contour is initialized by combining a convex hull shape with known parts (see Fig. 2b). Once the contour is deformed, it is combined with the extracortical borders found in the previous step and the final contour is saved. 2.2. Matching the extracted contour to the generic models A number of publicly available femur models were used to create the library of generic models: (i) the Standardized Femur [31], (ii) a model of a later generation composite femur (Pacific Research Labs, Vashion Island, WA) [32], (iii) models from the VAKHUM dataset—right and left femora [33], and (iv) a model from the Visible Human dataset [34]. According to a study by Noble et al. [28], human femur length varies between 353 and 557 mm (observed in a dataset of 200 patients) and correlates well with femoral head diameter and femoral neck length, which, in turn, correlate well with each other. The three models in the library were scaled to a reference length equal to the mean femur length from the data of Noble et al. (447.07 mm). The first step in matching the extracted contour to the generic models is to align the contour with each of the models in the same manner so as to identify the closest one. This is done by computing a Least Squares line of the medullary canal for generic models in the proximal part of the femur and similarly for the contour. Secondary axes are computed then by fitting Least Squares lines to the femoral neck points. Then these lines are merged so that the intersection point of the contour axes coincides with a projected intersection point of a 3D model. Taking into account cross-correlations of the femur dimensions [28], generic models are first scaled to the patient’s femur dimensions using the ratio between average femoral head diameter (46.59 mm) and length (447.07 mm). This step avoids the need to develop a method for dealing with cross-correlation between

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188

1183

Fig. 3. Similarity measure using correlation coefficient: (left) cc = 0.899, (middle) cc = 0.784, and (right) cc = 0.639.

radiograph scaling and object depth since it estimates coronal and sagittal femur dimensions based on anatomical relationships proposed by Noble et al. [28]. These models are then rotated about the medullary canal alignment axis from −15 to +15 degrees with a step of 1 degree to simulate femur rotations during radiography in an attempt to account for the unknown alignment of a femur when radiographed. A projection of the model to the coronal plane for each step creates a binary image which is further used to compute similarity with a binary image of the contour. Correlation coefficients computed for these images with respect to the contour image gives a similarity measure (Fig. 3) with the largest positive value identifying the closest model and rotation angle. 2.3. Warping the generic model The scaling procedure in the previous section provides a reasonably close approximation of the periosteal sizes of femur parts. However, according to Noble et al. [28], femoral neck shaft angle does not correlate well with periosteal dimensions. Also, Barratt et al. [10] observed that variation of femur shape has two principal modes lined to the following dimensions: length and neck angle. A warping algorithm proposed in this study consists of two steps: (i) warping femoral neck and head, and (ii) warping the femoral shaft.

obtained from the radiograph and the one generated during the search for the closest model from the library are finely sliced with a step of 0.01 mm (distal sections below the most distal position of the radiograph contour are not considered). The coronal plane passing through the femur shaft axis splits these segments in two parts (Fig. 5). Thus, the slicing results in four sets of segments—two for the 3D model and two for the radiograph-based contour (the medial end to the intersection point and the intersection point to the lateral end, see Fig. 5). Dividing by the lengths of the corresponding medial segments of the other model gives a scaling ratio (lateral scaling ratios are also computed in this way). Each pair of ratios, medial and lateral, defines two ellipse halves (see Fig. 5, inset). Warping is then performed as follows: for each point on the source surface, a radial ratio is computed by projecting the point onto the 2D shape previously described and, in turn, this ratio defines the target position of the point in the same radial direction. When a source point does not lie in a slicing plane, the two closest 2D shapes are linearly interpolated to compute required ratios. Finally, a layer at the proximal end of the warped shaft is smoothed by applying a similar filter as given by Eq. (2).

2.3.1. Femoral neck and head First, femoral neck angles are measured for the closest model and the extracted contour. An ellipsoid bounding femoral neck and head is intersected with the femur geometry (Fig. 4). This region is then rotated by the difference in femoral neck angles. The rotation is performed about a vector, defined as a cross product of the direction vectors of femoral neck and shaft axes. However, this transformation distorts the surface of the femur neck, which must then be smoothed by a filter (Eq. (2)):

⎛



⎜ ⎝1−

(, )= ⎜



1

exp

˘(1 − )



1 ˘(1 − )

−

 exp

(−((1 − )/2))  −

2



(−((1 − )/2)) 



∈

, 2



, +

,

∈

+

1− 2



1− ,1 2

(2)

where  ∈ [0,1] is the relative distance from the centre of the ellipsoid, ␯ is the parameter describing thickness of a layer to be smoothed, i.e. a part of the femoral neck, and ␥ is the smoothing behaviour parameter, empirically found to be best at a value of 0.2. 2.3.2. Femoral shaft Next, the femoral shaft up to the greater trochanter is deformed to match the patient’s geometry. For this purpose, contours

Fig. 4. Schematic explanation of bounding ellipsoid for femoral head and femoral neck: side view (left) and top view (right).

1184

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188

Fig. 5. Warping scheme for distal femur: slicing procedure in coronal planes defines medial and lateral segments; (inner) scaling radii define radial displacements med = gm gm c gm gm c c c , lat = Slat /Slat , where Smed , Smed , Slat , Slat are the lengths of medial and lateral segments for contour and generic model. Smed /Smed

2.4. Robustness test and geometric validation

2.5. Mechanical validation

Three tests were devised to test the reconstruction algorithm: a robustness test, a self-reconstruction test, and, finally, a comparison against the gold standard reconstruction method, 3D CT scans. First, 32 radiographs were selected from a dataset provided by the Adelaide & Meath Hospital, Tallaght, Dublin, Ireland. These radiographs were selected to contain proximal femur and pelvis, have no non-physiological objects, and to have no major contrast/brightness/position problems. All these radiographs were subjected to the reconstruction algorithm to assess robustness and applicability. Secondly, a model from the library of generic models was taken, a projection to the coronal plane was computed (virtual radiograph), and then this contour was used to find the closest generic model and warp it; this would allow measurement of error when the algorithm was used to reconstruct a model from its own projection. Next, the model used to create the virtual radiograph was removed to test how the algorithm performed when there is not a close match in the library of generic models. Reconstructed geometries were compared to the originals using a distance map approach (see Eq. (3)):

To understand how error in geometry may affect finite element analysis (FEA) of a proximal femur generated using the proposed reconstruction technique the following approach was adopted. A scheme proposed by Viceconti et al. [35] and later used by Ramos and Simões [36] was used to compare several FE models with varying degrees of error. First, a CT-reconstructed femur model was used as a reference as it is likely to contain the smallest error. Secondly, a model reconstructed from a corresponding radiograph was used as a model with introduced error. Finally, a model extrapolated from the previous two was created as the one containing largest error. These models were meshed using a mapping scheme so that the number of nodes, elements, and connectivity (and hence degrees of freedom) and their approximate positions were preserved. To investigate sensitivity to geometric inaccuracy, intermediate models were generated by linearly interpolating the original models. All models were restrained at the distal end and loads proposed by Viceconti et al. [35], representing joint and gluteus maximus forces (1976 and 1240 N, respectively), were applied at the femoral head and greater trochanter on the surface of the bone. A linear-elastic homogeneous material property was used for the femurs to remove the influence of variable material properties and hence focus on the influence of geometric error, similar to Viceconti et al. [35]; Young’s modulus and Poisson’s ratio were 14,200 MPa and 0.3, respectively. Linear-elastic analysis was performed using Abaqus (SIMULIA, Providence, RI, USA). For each model ten virtual gauges on medial and lateral sides of the femur were used to compare strain distributions.

d() ¯ = min dist(, ¯ T) ∀T ∈ St

(3)

where T is the single facet in a triangulated representation of the target surface, St , and ¯ is the point on the source surface. Finally, four patient datasets containing pelvic CT scans and radiographs for each patient were provided by the Adelaide & Meath Hospital, Tallaght, Dublin, Ireland. CT scans (512 × 512 pixels, 30–70 slices, 3 mm spacing) were used to create both right and left reference femur models using Slicer (Harvard Medical School, Boston, MA, USA). The planar radiographs for each patient (2136 × 1760 pixels, magnification factor of 5) were subjected to the proposed reconstruction algorithm. Radiograph-reconstructed models were aligned with the CT-reconstructed models using an Iterative Closest Point transform and were then compared to each other using the distance map approach described by Eq. (3).

3. Results 3.1. Robustness test and geometric validation Thirty-two images were successfully processed by the reconstruction method. An average running time was approximately 5 min on a computer with Intel Pentium® 4 CPU, 3 GHz processor. Approximately 90% of the runtime was consumed by the active contour algorithm. Results for reconstruction quality showed that the average distance from original to reconstructed model in case of the self-

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188 Table 1 Mean and standard deviation of distances between X-ray reconstructed models and CT reconstructed models for each patient of the CT-radiograph dataset. Patient ID

Right femur

Patient#1 Patient#2 Patient#3 Patient#4

3.283 2.573 3.743 3.101

± ± ± ±

Left femur

0.04 0.02 0.06 0.03

2.931 3.111 2.739 2.876

± ± ± ±

0.03 0.03 0.04 0.04

reconstruction was 0.48 ± 0.001 mm. For the case when the model to be reconstructed was excluded from the library of generic models (i.e. analogous to a reconstruction when no library model is a close match) this value increased to 1.44 ± 0.02 mm. For the four datasets with both CT scans and X-rays, reconstructed femora had an error in geometry between 2.57 and 3.74 mm (Table 1). Regions with the peak error (approximately 16 mm) were located in the femoral head (Fig. 6). However, if this problematic region was ignored for the calculation of the distance map, analogous to a hip replacement pre-operative plan, mean error value dropped down by approximately 0.5 mm (i.e. mean distance errors of between 2 and 3.25 mm). 3.2. Mechanical effect of geometric quality For the set of original and interpolated/extrapolated femora the maximum surface error ranged from 0.44 to 11.66 mm (Table 2). The mean surface error varied between 0.05 and 2.88 mm (Table 2). A geometric error coefficient was used to represent the level of interpolation/extrapolation ranging between 0 and 2, where 0 is the original model, 1 is the radiograph-based reconstructed, and 2 corresponds to the extrapolated model (Table 2 and Fig. 7). The largest surface error observed in extrapolated models stands in line with the error level observed in the models used to compare reconstruction quality from CT scans and planar radiographs. Measurement of octahedral shear strain at gauges on medial and lateral sides showed that error in octahedral shear strain was predicted to increase almost linearly with increasing geometric error (Fig. 7a). This near linearity enabled all further results to be summarised using three models: original, reconstructed and extrapolated. For the reconstructed model average strain errors for maximum, middle and minimum principal strain components

1185

were 30.228 ± 38.916, 13.593 ± 15.137, and 32.981 ± 46.013 ␮␧, respectively. For the extrapolated model these values were 55.917 ± 72.479, 24.982 ± 27.739, and 61.076 ± 84.332 ␮␧. Similar to the previous results, proportional growth in strain error was predicted with increase of geometric error. Examination of contour plots of strain error shows that the greatest error appears to lie on lateral and medial sides of the distal femur (see Fig. 7b).

4. Discussion Successful application of the proposed technique to 32 planar pelvic radiographs showed that it is robust, rapid and does not require user intervention, i.e. it is fully automated. Models reconstructed from corresponding CT scans and planar radiographs consistently contained a mean distance error of less than 3.7 mm. Peak errors were most often found to be in the regions of the femoral head and the greater trochanter. Additionally, error values doubled when there is not a relatively close match in the library of generic models. Investigation of the influence of geometric error on strain predictions revealed an almost linear relationship between predicted strain error and reconstructed geometry error for the range of error investigated. Finite element analysis performed on the models reconstructed using the proposed algorithm suggests that strain error caused by the reconstruction process can be relatively small. Results of the present study showed that the technique developed to automatically generate surface models from planar images was suitable to perform the task. Evaluation of accuracy was carried out and showed that average surface distance error was equal to 3.04 mm; this result stands in line with error magnitudes computed using other reconstruction techniques, which were within the following range—0.7–6.88 mm [10–12, 16, 19, 21]. One benefit of the technique proposed in this study is that it is highly automated and able to generate surface models rapidly. This can be used to process large amounts of patient-specific radiographs for morphological analysis, measuring anatomical variations in proximal femur shape as well as pre-operative planning simulations. Investigation of the influence of geometric error on the mechanical response of the loaded bone showed that growth of the geometric error results in a proportional increase in strain error magnitude. It also indicated that change in precision up to the level of generated models may not change the strain picture dramatically (Fig. 7). Therefore,

Fig. 6. Example error distributions for three patients from the CT–X-ray reconstruction comparison tests (numbers stated for each patient are the mean error). Peak surface errors occurred at the proximal surface of the femoral head in all models. Table 2 Surface error between tested and reference models, mm: dmax is the maximum observed error, dmean is the average distance, and EC is the linear coefficient representing error level. Model ID

[1]1

2

3

4

5

6

7

8

dmax dmean EC

0.437 0.048 0

1.456 0.370 0.25

2.913 0.727 0.5

4.370 1.085 0.75

5.827 1.443 1

7.284 1.801 1.25

8.741 2.159 1.5

10.198 2.517 1.75

9 11.656 2.876 2

1186

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188

Fig. 7. (a) Schematic positions of strain gauges and strain measurements for different surface error levels. Error coefficient indicates degree of error in model: 0, CT-based reconstruction; 1, radiograph-based reconstruction; and 2, extrapolated model based on parameters estimated from linear interpolation between model 0 and 1. (b) Strain error contour plots (relative to the CT-reconstructed model).

models reconstructed using the proposed 2D–3D reconstruction technique may be used to carry out patient-specific analysis with an estimate of the likely error when CT geometry is not available. Possible limitations of the proposed method are low-quality radiographs (e.g. if an uncalibrated X-ray machine was used to acquire images), artefacts generated by presence of metallic and other objects that could mislead the contour recognition tech-

nique, and femur shapes which are not covered by the assumed size distribution, e.g. extreme size of different ethnicities or anatomic deviations such as high grade dysplastic femora. Also, correlation analysis of the generic model projection contour and the X-ray extracted contour was used to select the best matching model from the library of generic models. More sophisticated procedures, such as use of Fourier descriptors as proposed by Messmer et al. [37], may

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188

select a better match that could improve the subsequent operations to give a more accurate final model; further investigation comparing the techniques could be undertaken in future to assess if this is worthwhile. Another limitations is that planar pelvic radiographs do not capture the three-dimensional morphology of the greater trochanter region very well, in particular the attachment location of the piriformis muscle. Thus, this technique could be improved for reconstruction of this region. Although an attempt was made to account for the unknown orientation of the femur in the radiograph using an iterative rotation–projection–correlation procedure, no sensitivity analysis was performed to assess what the impact of unknown orientation on the reconstruction technique might be. One possible method of testing this would be to rotate a model or excised femur over a range of angles, radiographing it for each angle, and then reconstructing it and comparing the reconstructions with the parent femur. Finally, for the current investigation, only the most proximal part of the femur (from just below the lesser trochanter up) was used to validate the accuracy against the CT scans—it is expected that this region contains the largest errors from the reconstruction procedure, thus, making the erroneous part proportionally large compared to other investigations. From the surface error distributions of the CT–X-ray pairs, one can see that the femoral head was a source of inaccuracy for reconstructed geometries. If the femoral head was removed and distances were measured again, the average error magnitude decreased from 3.04 to 2.31 mm. Analyzing visual results in Zheng et al. [19] and Barratt et al. [10] it is possible to notice similar problems. This may happen due to various reasons, one of which is an insufficient resolution of either the CT image or the X-ray image. Analysis of the CT slices containing the femoral head showed that the vertical spacing is too large to reconstruct a precise geometry. In spite of this limitation, surfaces reconstructed using the proposed algorithm may still be useful for applications such as pre-operative patient-specific simulations of total hip replacements, where the femoral head is being removed and replaced by an implant. For the FEA validation, predicted strain and geometric error appear to be almost linearly proportional. However, there could be several limitations of the study that caused such a prediction of near linear change in strain error with changing geometric error. One possibility is the use of a common mesh mapped to the different geometries. This was done to maintain the same mesh topology, virtual strain gauge placement, and number of degrees of freedom across all models. However, a side effect of mapping the mesh in this way is that erroneous regions that grow in size have larger elements than corresponding regions in the baseline model and thus reduced accuracy in those regions. This could have the effect of smearing out a nonlinear change in strain predicted in such a region, causing it to be perceived as linear. Another factor could be placement of the gauges, most of which were placed along the medial and lateral sides of the femur meshes. This was chosen due to the loading, which placed the femur primarily in bending, so that the medial and lateral shaft surfaces would be expected to have the highest sensitivity to geometric error. However, these surfaces tended to have the lowest errors since medial and lateral borders were extracted directly, rather than inferred, from the radiographs and the femur shape is simplest along the shaft. Thus, the geometric error range was relatively low in these regions so that even a relatively nonlinear response could appear linear because of the small field of view of the response. More accurate application of physiological loads may be required, e.g. more distributed muscle loading (a region instead of selected nodes) and a larger number of muscles and activities (to make regions with higher geometric error more sensitive to strain error), to improve understanding of how the geometric accuracy affects strain prediction in femur models and to determine whether the result of linearity might be generalised to other cases.

1187

In spite of these limitations, this study shows clearly that a patient-specific surface model can be reconstructed automatically from a planar radiograph image when volumetric data is not available, and the reconstructed geometry can be used in a number of orthopaedic applications. Acknowledgements This study was funded by a grant from Science Foundation Ireland under its Research Frontiers Programme (05/RFP/ENG0097). The authors would also like to thank Mr. Ruairi MacNiochaill, Dept. Orthopaedic Surgery, Cappagh National Orthopaedic Hospital, Ireland and the radiology department of the Adelaide and Meath Hospital, Tallaght, Dublin, Ireland for help in obtaining CT scans and X-rays used in this study. Conflict of interest statement: There is no conflict of interest in relation to this work. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.medengphy.2010.08.009. References [1] McCollough CH, Primak AN, Braun N, Kofler J, Yu L, Christner J. Strategies for reducing radiation dose in CT. Radiologic Clinics of North America 2009;47(1):27–40. [2] Krüger J, Westermann R. Acceleration techniques for GPU-based volume rendering. In: Tork G, van Wijk JJ, Moorhead R, editors. VIS 2003. Seattle Washington, USA: IEEE; 2003. p. 287–92. ISBN: 0-7803-8120-3, doi:10.1109/VISUAL.2003.1250384. [3] Keyak JH, Meagher JM, Skinner HB, Mote Jr CD. Automated three-dimensional finite element modelling of bone: a new method. Journal of Biomedical Engineering 1990;12(5):389–97. [4] Lorensen WE, Cline HE. Marching cubes: a high resolution 3D surface construction algorithm. Computer Graphics 1987;21(4). [5] Hounsfield GN. Nobel lecture, 8 December 1979. Computed medical imaging. Journal of Radiology 1980;61(6–7):459–68. [6] Couteau B, Hobatho M-C, Darmana R, Brignola J-C, Arlaud J-Y. Finite element modelling of the vibrational behaviour of the human femur using CT-based individualized geometrical and material properties. Journal of Biomechanics 1998;31(4):383–6. [7] Yosibash Z, Trabelsi N, Milgrom C. Reliable simulations of the human proximal femur by high-order finite element analysis validated by experimental observations. Journal of Biomechanics 2007;40(16):3688–99. [8] Viceconti M, Davinelli M, Taddei F, Cappello A. Automatic generation of accurate subject-specific bone finite element models to be used in clinical studies. Journal of Biomechanics 2004;37(10):1597–605. [9] Viceconti M, Zannoni C, Testi D, Cappello A. CT data sets surface extraction for biomechanical modeling of long bones. Computer Methods and Programs in Biomedicine 1999;59(3):159–66. [10] Barratt DC, Chan CSK, Edwards PJ, Penney GP, Slomczykowski M, Carter TJ, et al. Instantiation and registration of statistical shape models of the femur and pelvis using 3D ultrasound imaging. Medical Image Analysis 2008;12(3):358–74. [11] Rajamani KT, Styner MA, Talib H, Zheng G, Nolte LP, Ballester MAG. Statistical deformable bone models for robust 3D surface extrapolation from sparse data. Medical Image Analysis 2007;11(2):99–109. [12] Talib H, Rajamani K, Kowal J, Nolte LP, Styner M, Ballester MAG. A comparison study assessing the feasibility of ultrasound-initialized deformable bone models. Computer Aided Surgery 2005;10(5–6):293–9. [13] Heimann T, Meinzer H-P. Statistical shape models for 3D medical image segmentation: a review. Medical Image Analysis 2009;13(4):543–63. [14] Fleute M, Lavallee S, Julliard R. Incorporating a statistically based shape model into a system for computer-assisted anterior cruciate ligament surgery. Medical Image Analysis 1999;3(3):209–22. [15] Shim VB, Pitto RP, Streicher RM, Hunter PJ, Anderson IA. The use of sparse CT datasets for auto-generating accurate FE models of the femur and pelvis. Journal of Biomechanics 2007;40(1):26–35. [16] Bras AL, Laporte S, Bousson V, Mitton D, Guise JAD, Laredo JD, et al. Personalised 3D reconstruction of proximal femur from low-dose digital biplanar radiographs. International Congress Series 2003;1256:214–9. [17] Zanetti EM, Crupi V, Bignardi C, Calderale PM. Radiograph-based femur morphing method. Medical & Biological Engineering & Computing 2005;43(2):181–8. [18] Zheng G. Statistically deformable 2D/3D registration for accurate determination of post-operative cup orientation from single standard X-ray radiograph.

1188

[19]

[20]

[21]

[22]

[23]

[24]

[25] [26] [27]

P.E. Galibarov et al. / Medical Engineering & Physics 32 (2010) 1180–1188

Medical Image Computing and Computer-Assisted Intervention 2009;12(Pt 1):820–7. Zheng G, Gollmer S, Schumann S, Dong X, Feilkas T, Ballester MAG. A 2D/3D correspondence building method for reconstruction of a patient-specific 3D bone surface model using point distribution models and calibrated X-ray images. Medical Image Analysis 2009;13(6):883–99. Zheng GY. Statistical shape model-based reconstruction of a scaled, patientspecific surface model of the pelvis from a single standard AP X-ray radiograph. Medical Physics 2010;37(4):1424–39. Langton CM, Pisharody S, Keyak JH. Comparison of 3D finite element analysis derived stiffness and BMD to determine the failure load of the excised proximal femur. Medical Engineering & Physics 2009;31(6):668–72. Langton CM, Pisharody S, Keyak JH. Generation of a 3D proximal femur shape from a single projection 2D radiographic image. Osteoporosis International 2009;20(3):455–61. Bredbenner TL, Nicolella DP. Reconstructing High Fidelity 3D Finite Element Models of a Human Proximal Femur from 2D DXA Data. In: Transactions of the 54th Annual Meeting of the Orthopaedic Research Society. The Orthopaedic Research Society; 2008. p. Poster No. 911. Lennon AB, Britton JR, MacNiocaill RF, Byrne DP, Kenny PJ, Prendergast PJ. Predicting revision risk for aseptic loosening of femoral components in total hip arthroplasty in individual patients—a finite element study. Journal of Orthopaedic Research 2007;25(6):779–88. Hou Z, Han C. Force field analysis snake: an improved parametric active contour model. Pattern Recognition Letters 2005;26(5):513–26. Caselles V, Kimmel R, Sapiro G. In: Bovik A, editor. Geometric active contours for image segmentation. Burlington: Academic Press; 2005. p. 613–27. Lim JS. Two-dimensional signal and image processing. Prentice Hall: Englewood Cliffs, NJ; 1990.

[28] Noble PC, Alexander JW, Lindahl LJ, Yew DT, Granberry WM, Tullos HS. The anatomic basis of femoral component design. Clinical Orthopaedics and Related Research 1988;(235):148–65. [29] Sugano N, Noble PC, Kamaric E. Predicting the position of the femoral head center. The Journal of Arthroplasty 1999;14(1):102–7. [30] Leavers VF. Which Hough transform? CVGIP: Image Understanding 1993;58(2):250–64. [31] Viceconti M, Casali M, Massari B, Cristofolini L, Bassini S, Toni A. The ‘Standardized femur program’ proposal for a reference geometry to be used for the creation of finite element models of the femur. Journal of Biomechanics 1996;29(9):1241–1241. [32] Cheung G, Zalzal P, Bhandari M, Spelt JK, Papini M. Finite element analysis of a femoral retrograde intramedullary nail subject to gait loading. Medical Engineering & Physics 2004;26(2):93– 108. [33] Van Sint Jan S. The VAKHUM—project: virtual animation of the kinematics of the human. Theoretical Issues in Ergonomics Science 2005;6(3):277– 1277. [34] Ackerman MJ. The Visible Human Project. The Journal of biocommunication 1991;18(2):14–14. [35] Viceconti M, Bellingeri L, Cristofolini L, Toni A. A comparative study on different methods of automatic mesh generation of human femurs. Medical Engineering & Physics 1998;20(1):1–10. [36] Ramos A, Simões JA. Tetrahedral versus hexahedral finite elements in numerical modelling of the proximal femur. Medical Engineering & Physics 2006;28(9):916–24. [37] Messmer P, Long G, Suhm N, Regazzoni P, Jacob AL. Volumetric model determination of the Tibia based on 2D Radiographs using a 2D/3D database. Computer Aided Surgery 2001;6(4):183–94.