Computer Vision and Image Understanding 113 (2009) 202–211

Contents lists available at ScienceDirect

Computer Vision and Image Understanding
journal homepage: www.elsevier.com/locate/cviu

3D reconstruction of a femoral shape using a parametric modeland two 2D ﬂuoroscopic images
Ryo Kurazume a,*, Kaori Nakamura a, Toshiyuki Okada b, Yoshinobu Sato c, Nobuhiko Sugano c,
Tsuyoshi Koyama d, Yumi Iwashita a, Tsutomu Hasegawa aaGraduate School of Information Science and Electrical Engineering, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan
Osaka University Hospital, Medical Center for Translational Research, 2-15 Yamadaoka, Suita-shi, Osaka 565-0871, Japanc
Graduate School of Medicine, Osaka University, 2–2 Yamadaoka, Suita-shi, Osaka, 565-0871, Japand
Department of Orthopaedic Surgery, National Hospital Organization Osaka-Minami Medical Center, 2-1 Kidohigashimachi, Kawachinagano, Osaka 586-8521, Japanba r t i c l e

i n f oArticle history:
Received 18 September 2007
Accepted 15 August 2008
Available online 12 September 2008
Keywords:
Fluoroscopic image
Parametric femoral model
Registration
Medical image diagnosisa b s t r a c t
In medical diagnostic imaging, the X-ray CT scanner and the MRI system have been widely used to examine 3D shapes and internal structures of living organisms and bones. However, these apparatuses are generally large and very expensive. Since an appointment is also required before examination, these systems
are not suitable for urgent fracture diagnosis in emergency treatment. However, X-ray/ﬂuoroscopy has
been widely used as traditional medical diagnosis. Therefore, the realization of the reconstruction of precise 3D shapes of living organisms or bones from a few conventional 2D ﬂuoroscopic images might bevery useful in practice, in terms of cost, labor, and radiation exposure. The present paper proposes amethod by which to estimate a patient-speciﬁc 3D shape of a femur from only two ﬂuoroscopic images
using a parametric femoral model. First, we develop a parametric femoral model by the statistical analysis of 3D femoral shapes created from CT images of 56 patients. Then, the position and shape parametersof the parametric model are estimated from two 2D ﬂuoroscopic images using a distance map constructed by the Level Set Method. Experiments using synthesized images, ﬂuoroscopic images of a phantom femur, and in vivo images for hip prosthesis patients are successfully carried out, and it is veriﬁedthat the proposed system has practical applications.
Ó 2008 Elsevier Inc. All rights reserved.

1. Introduction
In medical diagnostic imaging, the X-ray Computed Tomography (CT) scanner and the Magnetic Resonance Imaging (MRI) system have been widely used to examine the 3D shape or internal
structure of living organisms and bones. However, these apparatuses are generally large and very expensive, and thus, they areusually installed in large medical institutions rather than small
local clinics. Since an appointment is also required before examination, these systems are not suitable for urgent fracture diagnosis inemergency treatment.
Meanwhile, X-ray/ﬂuoroscopy has been widely used as traditional medical diagnosis. Recently, digital ﬂuoroscopy has beendeveloped and widely used in many hospitals. The cost of this ﬂuoroscopic inspection system is much lower than that of CT or MRI
systems and the system can be dealt with more conveniently. Furthermore, the risk of radiation exposure is also lower than that of
the CT inspection system.
From the above considerations, the realization of the reconstruction of precise 3D shapes of living organisms or bones from a few* Corresponding author. Fax: +81 92 802 3607.
E-mail address: kurazume@is.kyushu-u.ac.jp (R. Kurazume).
1077-3142/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved.
doi:10.1016/j.cviu.2008.08.012

conventional 2D ﬂuoroscopic images might be very useful in practice, in terms of cost, labor, and radiation exposure. In particular,
there is a strong demand from surgeons for 3D computer-aided surgery without laborious CT imaging for simple surgeries such as artiﬁcial joint replacement or fracture treatment. A practical 3D
diagnostic system using common 2D ﬂuoroscopic images is desired.
However, 3D shape reconstruction from a 2D image is a fundamentally ill-posed problem, and so a sufﬁcient number of images
must be obtained, or several constraint conditions for the 3D shapemust be determined. However, the shapes of bones have inherent
and universal patterns, and thus, by modelling such inherent patterns, 3D shape reconstruction from a few 2D images is possible.
In the present paper, a technique by which to estimate thepatient-speciﬁc 3D shape of a femur from only two ﬂuoroscopic
images is proposed. The proposed technique utilizes a parametricfemoral model constructed by statistical analysis of 3D femoral
shapes created from CT images of 56 patients. The position/orientation and shape parameters of the parametric model are then estimated from two 2D ﬂuoroscopic images by solving the 2D/3D
registration problem using a distance map constructed by the Level
Set Method.
The 2D/3D registration problem is well established in image
processing, especially for texture mapping in Computer Graphics

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211

or Augmented Reality. For a rigid object, (1) feature-based
techniques [1–3], (2) image-based techniques using texture, reﬂectance, brightness, and shading [4,5,6], and (3) silhouette-based
techniques [7–10], have been proposed. In particular, in surgicalnavigation systems, Digitally Reconstructed Radiographs (DRRs)[11,12] are widely used in 2D/3D registration for the ﬂuoroscopy-guided surgery. However, the construction of DRR is timeconsuming and special techniques such as the use of graphics
hardwares is indispensable for quick medical diagnosis. Since texture or shadow are usually unavailable in ﬂuoroscopic images, fast
2D/3D registration technique using simple and robust features
such as silhouette and contour lines is desirable.
In 2D/3D registration of a non-rigid object such as soft tissues inmedical imaging, similarity measure [13,14], mutual information[15], afﬁne [16,17], geometric hashing [18], and displacementﬁeld-based transformation [19] have been proposed and tested.
In addition, the 3D shape estimation of a parameterized object,such as the shape reconstruction of mathematical plaster models
with unknown parameters using a laser range ﬁnder [20], or thecomparison of multiple cross-section images of a 3D model and a
3D parametric model [21], has also been studied. However, thesestudies assumed the use of a sufﬁcient number of images or a precise 3D shape taken by a laser range ﬁnder, and only a few studies
have examined 3D non-rigid shape reconstruction from only a few
2D images [22–24]. Zheng et al. [24] proposed a similar approachwith our method for estimating a femoral shape from ﬂuoroscopic
images. Their technique is based on an active shape model andconventional ICP method. They introduced experiments for 11
cadaveric femurs using three ﬂuoroscopic images including theone taken from the longitudinal direction of the femoral shaft.
However, detailed discussion using in vivo images taken by clinicalsetting and calibration technique for ﬂuoroscopic imaging was notpresented.
This paper propose the 3D reconstruction technique of a femoral shape using two in vivo ﬂuoroscopic images of patients. A distance map constructed by the Level Set Method is utilized for fastand robust 2D/3D registration with M-estimator. The calibration
technique for ﬂuoroscopic imaging using two kinds of calibration
markers is also presented. This paper is organized as follows. Section 2 describes the statistical shape model of the femur proposedby Okada, et al. [25,26] at ﬁrst. Then, we introduces the proposed
3D reconstruction technique of a femoral shape using the statistical shape model and two 2D ﬂuoroscopic images. In Section 3,
experiments using synthesized images, ﬂuoroscopic images of aphantom femur, and in vivo images for hip prosthesis patients,
are successfully carried out, and it is veriﬁed that the proposed system has practical applications. Finally, Section 4 presents ourconclusion.

2. 3D parametric femoral model
In this section, we introduce the 3D reconstruction technique ofa femoral shape using the statistical shape model and two 2D ﬂuoroscopic images.
2.1. Overview of the proposed algorithm
The procedure of the proposed technique is summarized asfollows:(1) First, we develop a 3D parametric femoral model by the statistical analysis of 3D femoral shapes created from CT
images of 56 patients. With this model, a general 3D shapeof the femur is expressed by the average shape and several
shape parameters. Thus, the 3D shape estimation of the203

patient’s femur from two 2D ﬂuoroscopic images can bedivided into two procedures, the determination of optimum
position and orientation of the 3D model in two 2D ﬂuoroscopic images (Step 2) and the estimation of optimum shapeparameters (Step 3).
(2) The position and orientation of the parametric model in two
2D ﬂuoroscopic images are determined by the 2D/3D registration technique using a distance map constructed by theLevel Set Method.
(3) The optimum shape parameters of the parametric model are
estimated by comparing the silhouette contour of the parametric model and two 2D ﬂuoroscopic images using a distance map obtained in Step 2.
(4) Step 2 and 3 are repeated until the residual error between
the silhouette contour of the parametric shape model andtwo 2D ﬂuoroscopic images becomes less than a thresholdvalue.
Each of the above steps is explained in details in the following
sections.
2.2. 3D parametric femoral model
We utilize the statistical shape model of the femur proposed by
Okada et al. [25,26]. In this technique, a number of 3D femoral
shapes created from CT images are analyzed statistically, and theparametric femoral model [27,28], which consists of the average
3D shape and several shape parameters, is created. With this parametric femoral model, a general 3D shape of the femur is expressedby the average shape and several shape parameters.
The concrete procedure for creating a parametric 3D femoral
model is as follows:
1. Surface models of femurs are created from CT images by manual segmentation and Marching Cubes.
2. Local coordinate axes of the surface models are determined byapplying the principal component analysis (PCA) to the set of
3D positions of the node points in each surface model. The center of gravity is deﬁned as the origin of the local coordinate system. The Z axis is determined as the axis corresponding to thelargest eigenvalue, which is toward the longitudinal directionof the femoral shaft from the hip to knee. The region up to
35% of total length of the femur along the Z axis from the femoral head top deﬁned as the node point having the minimum Z
coordinate is extracted as a proximal femur, where the total
length is deﬁned as the difference of the maximum and minimum Z coordinates in the femur surface model, and the valueof 35% is determined so as to cover a region anatomicallyregarded as a proximal femur.
3. One of the femoral model is selected as the reference model anddisplacement vector ﬁelds to all other models described by thethin-plate spline are calculated using a point-based non-rigid
registration algorithm[29]. Although the approximation accuracy of the resulted parametric model may somewhat dependon the selected reference model, we currently selected it basedon visual assessment so that it is not largely deviated from theaverage shape.The reference model is decimated using an algorithm available in the visualization toolkit so that the number ofnodes of the surface model becomes 1500. 1500 points wereused because the proximal femur shape can be represented ina reasonable accuracy using 1500 points while computational
cost for the non-rigid registration is still acceptable, for example, less than 10 min for each case.The reference model is
non-rigidly registered with each of other models so as to ﬁnd
the corresponding 1500 points on each model. The non-rigid
registration algorithm [29] generate 3D displacement vector

204

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211

ﬁeld from the reference to each dataset. The 1500 points of eachmodel are determined by combining those of the reference
model and the generated displacement vector ﬁeld. These processes do not guarantee physically-meaningful correspondences. Using these correspondences, however, each surface
model can be represented as a 1500  3 dimensional shape vector, that is, a ﬁxed-dimensional vector. For our current purpose,correspondences are acceptable if they are plausible eventhough there is not some guarantee.
4. Given n shape vectors, the average shape vector is given by theiraverage. PCA is applied to a set of the shape vectors subtracted
by the average shape vector to obtain the eigenvectors whose
coefﬁcients correspond to the shape parameters.
The parametric femoral model used in the following experiments was created using CT images of 56 patients. By applying
PCA to 56 samples of 3D femoral shapes, we extracted the mostsigniﬁcant 50 principal components (p1 ; p2 ; . . . ; p50 ), standard deviation (r1 ; r2 ; . . . ; r50 ), and corresponding principal vectors(v1 ; v2 ; . . . ; v50 ). With the obtained parametric femoral model, thegeneral 3D shape of a femur is expressed asx0 ¼ x þ ðp1  r1  v1 Þ þ ðp2  r2  v2 Þ þ   

ð1Þ
0where x is the surface point of the average shape and x is the surface point of the general shape. Therefore, the general 3D shape of afemur is expressed by the parametric femoral model with
 average 3D shape and several principal vectors
determined) several (up to 50) shape parameters (estimated)(pre-Fig. 1 shows the contribution ratio of the shape parameters forthe statistical femoral model.
2.3. Reconstruction of 3D femoral shape from two 2D ﬂuoroscopic
imagesIn this section, we introduce the 2D/3D registration algorithm
and the estimation procedure of the optimum shape parameters
using two ﬂuoroscopic images.
This 2D/3D registration algorithm utilizes the contour lines ofthe silhouette of the 2D image and the projected contour lines of
the 3D model. The optimum position of the 3D model is determined such that the contour lines coincide with each other on
the 2D image plane. In commonly used approaches such as the
ICP algorithm, the error metric is usually deﬁned as the sum of193%

the distances between the points on the 2D contour lines and theirnearest points on the projected contour lines of the 3D model.
However, the nearest point search is a laborious task and is time
consuming even for the kd tree-based algorithm [30].
In the present approach, the 2D distance map [7] is utilized.
First, the 2D distance map from the contour lines is created on
the 2D image using the Fast Marching Method [31,32] or rasterscan algorithms using local operators [33,34]. Once the 2D distance
map is created, the error metric is obtained directly from the 2D
distance map as the value at the points on the projected contour
lines of the 3D model. Using the course-to-ﬁne strategy called‘‘Distance Band” [7], a 2D distance map can be constructed quite
rapidly using the Fast Marching Method.
When 2D/3D registration and estimation of the shape parameters are performed at the same time, the depth from the view point
and the scale of the 3D model cannot be distinguished. Therefore,
the proposed algorithm utilizes two ﬂuoroscopic images takenfrom two viewpoints at different positions. In addition, we assume
that the 3D femoral parametric model is constituted by a large
number of small triangle patches of approximately the same size.
In the following sections, we explain Step 2 and Step 3 shown in
Section 2.1. In both steps, 2D distance maps from the extracted
contour lines of patient’s femur in ﬂuoroscopic images are created
at ﬁrst. Then the error between the contour lines of the femur andthe projected contour lines of the 3D model is minimized by reading the distance values of the 2D distance maps.
2.3.1. Registration of 2D ﬂuoroscopic images and the 3D parametricmodel
A brief description of the registration procedure of the 2D ﬂuoroscopic images and the 3D parametric model is given as follows:
1. Extract the contour lines of the femur in the ﬂuoroscopic
images using an active contour model such as snakes or the
Level Set Method [31].
2. Construct a 2D distance map from the extracted contour lines
using the Fast Marching Method [31,32]. Fig. 2 shows an example of the constructed 2D distance map of a femoral image.
3. Place the parametric femoral model at an arbitrary initial position which is determined manually and calculate the 2D projection image of the 3D model.
4. Extract contour lines of the projected image and corresponding
3D patches of the 3D model. This procedure can be executed
quite rapidly by the OpenGL hardware accelerator (Appendix
A.1).
5(a). Apply the force calculated from the 2D distance map at the
projected contour points directly to the corresponding 3D
patch. Details are presented in Section 2.4.
6(a). Using the robust M-estimator, which is a robust estimation
technique, the total force and moment around the centerof gravity (COG) is calculated.

Contribution

0.8

0.6

+20 pixels+40 pixels

0.4

0.2

Boundary

0
0

10

20

30

Number of parameters
Fig. 1. Contribution of parametric model.

40

50

Contour line in distance map
Fig. 2. 2D distance map from contour in the femoral image.

205

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211

7(a). Steps 3 to 6(a) of the procedure are repeated for two images
captured from the different viewpoints sequentially, and thetotal force and moment are calculated.
8(a). Update the position of the 3D parametric model according to
the total force and moment.
9(a). Repeat Steps 3 to 8(a) until the magnitude of the total force
and moment becomes less than the pre-deﬁned thresholdvalue.
2.3.2. Estimation of the optimum shape parameters
The estimation procedure of the optimum shape parameters of
the 3D parametric femoral model is shown in this section. This procedure also uses the 2D distance map from the contour line of thefemur in the ﬂuoroscopic image, which has already been constructed, as described in Section 2.3.1. Therefore, from Steps 1 to
4 are same as Section 2.3.1 and we can skip these procedures byusing the obtained 2D distance map in Step 4 of the above procedure, After Step 4, the optimum shape parameters are estimated asfollows:
5(b). Calculate the error E, which is deﬁned as the sum of the values of the 2D distance map at the projected contour line of
the 3D parametric model.
6(b). Find the optimum shape parameters that minimize the error
E at the current position using the conjugate gradientmethod.
7(b). Reconstruct the 3D shape according the obtained shapeparameters using Eq. (1).
8(b). Repeat Steps 3 to 7(b) until the error E becomes less than the
pre-deﬁned threshold value.
2.4. 2D/3D registration using the robust M-estimator

ð2Þ

where Di is the value of the distance map at the contour point,
which corresponds to the triangular patch i, and rDi is the gradientof Di .
In Step 6(a), the total force and moment around the center ofgravity (COG) is calculated by the following equations:

F¼

X

M¼

wðf i Þ

i
X

ð3Þ

b
Center of gravity

ð4ÞCenter of gravity

riF
Objectfi

Mfi

Image plane
Projected
contour pointyxz

Steepest descentdirectionv
Focal point

Fig. 4. Calculation of the total force and moment around the center of gravity(COG).

where r i is a vector from the COG to triangular patch i and wðzÞ is aparticular estimate function. In practical situations, the contour ofthe femur is occasionally occluded or blurred, or the 2D image iscorrupted by noise. In such cases, the obtained boundary does notcoincide with the projected contour of the 3D model and the correct
distance value cannot be obtained. To deal with this problem, we
introduce the robust M-estimator in order to disregard contour
points with large errors. In our implementation, we utilized a following function using a Lorentzian function with a variance r2
and a vector z.

wðzÞ ¼

1
1þ j zj2 =r2zð5ÞLet us consider the force f i and the moment r i  f i as an error zi
and the sum of the error asXqðzi Þ

ð6Þ

iwhere EðPÞ is a scalar function of a vector P which is the position of
the 3D parametric model, and qðzÞ is a particular estimate function,
which is deﬁned asoqðzÞ
¼ wðzÞ
oz

ð7Þ

The position P that minimizes the error EðPÞ is obtained as thefollowing equation:oE X oqðzi Þ ozi¼¼0oPozi oPið8ÞHere, we deﬁne the weight function wðzÞ as the following equationin order to evaluate the error term:wðzÞz ¼ wðzÞ ¼

wðr i  f i ÞfifiEðPÞ ¼

After obtaining the distance map on the 2D ﬂuoroscopic image
and the list of the triangular patches of the 3D model corresponding to the contour points, the force f i is applied to all of the triangular patches of the contour points (Figs. 3 and 4), as explained in
Step 5(a).

rDi
f i ¼ Di
jrDi j

a

oqðzÞ
oz

ð9Þ

From the above equation, we obtain the following weighted leastsquares method:

i

oE X
ozi
wðzi Þzi
¼
¼0
oP
oP
i

Boundaryfi

ð10ÞIn our implementation, the optimum position, which minimizes theerror EðPÞ, is obtained by the steepest gradient method, as shown in
Step 9(a).
3. Experiments
3.1. Simulation using DRRs

Fig. 3. Force fi is applied to 3D triangular patch i on the contour line.

First, we conducted the experiments using Digitally Reconstructed Radiographs (DRRs) in order to evaluate the fundamental

206

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211

Ap

Fig. 5. Directions of DRRs for simulation.

repeated alternately and independently. An example of the experimental results for test data 4 is shown in Fig. 7, which illustrates
the average shape, the actual shape, and the estimated shape.
Fig. 8 indicates the average error between the estimated shapeand the actual shape. The average error is deﬁned as the average ofthe minimum distance from the surface point of the estimatedshape to the triangle patches of the actual shape. In this ﬁgure,‘‘AS” on the horizontal axis indicates the average error when theposition is estimated but all of the shape parameters are ﬁxed totheir initial values (all of the parameters are set to ‘‘0”). This ﬁgureshows that the average error gradually decreases as the number of
estimated shape parameters increases. However, the average error
converges when the number of shape parameters is approximately
ﬁve and no signiﬁcant difference is observed, even if the number ofthe shape parameters increases.
In addition, Fig. 9 shows the types of errors deﬁned below when
the number of estimated shape parameters is ﬁve. Table 1 also
indicates the average of the error, the standard deviation, the maximum value, and the minimum value for the 10 models.
Average error 1 The average error between the average modeland the actual shape.
Average error 2 The average error between the 3D optimized
estimated shape and the 3D actual shape obtained by comparing the actual shape and the parametric model.
Average error 3 The average error between the estimated shapeand the actual shape obtained by comparing the two 2D ﬂuoroscopic images and the 3D parametric model (proposed method).

performance of the proposed method. In the experiment, the estimation accuracy for 10 femoral models is examined using two
reconstructed ﬂuoroscopic images. Among 10 models, ﬁve models(model data 1–5) are used for the construction of the 3D parametric model and ﬁve models (test data 1–5) are not used.
We determine the directions of the ﬂuoroscopic images asshown in Fig. 5 considering the possible direction in actual radiographic examination. Under this condition, two view directions
meet at right angles at the main axis of the femur AP . Examplesof the reconstructed ﬂuoroscopic images are shown in Fig. 6.
All of the 3D femoral shapes used in the experiments werereconstructed precisely by the CT scanner beforehand, and the
optimized shape parameters, which minimize the distance errors
between surface points, were determined by comparing the 3D actual shape and the 3D parametric model and searching all possible
candidates.
First, we chose up to 10 principal components and estimatedthe position and optimum shape parameters of the femur on the
ﬂuoroscopic images. In this experiment, the position estimationof the femur and the optimum parameter estimation were

The experimental results show that Average error 3 between
the estimated shape and the actual shape is less than 1.1 mm at
2.5

test data 1~5
model data 1~5

2

Average error [mm]Fig. 6. Examples of DRRs for simulation experiments.

1.5

1

0.5

AS: Average shape
0

AS

2

4

6

10

Number of estimated parameters
Fig. 8. Average error and number of estimated shape parameters.

No.1

No.1
3

3

2
1
0

No.5

No.4

test data

No.2

No.3

average error 1
Fig. 7. Actual and estimated femoral model with and without shape parameter
estimation.

8

No.5

2
1
0

No.4

No.2

No.3
model data

average error 2

average error 3

Fig. 9. Average error after shape parameter estimation (number of principal
components: 5).

207

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211
Table 1Comparison of average errors (mm)Average error 1
Average error 2
Average error 3(Test data)(Model data)

Average

STD.

Maximum

Minimum

1.69
0.90

0.54
0.13

2.52
1.06

0.90
0.60

0.91
0.81

0.15
0.07

1.08
0.90

0.66
0.71

worst, and it is veriﬁed that the 3D shape can be estimated using
two 2D ﬂuoroscopic images with the same accuracy in case thatthe 3D shapes are compared directly. Moreover, it is conﬁrmedthat there is no signiﬁcant difference between the models thatwere used for the construction of the parametric femoral modeland those that were not.
Note that in some cases in Fig. 9, Average error 3 is smaller than
Average error 2. In the calculation of Average error 2, the origin andthe direction of the local coordinate systems of the estimated andactual shapes coincide with each other and relative position ofthese shapes is ﬁxed. On the other hand, for Average error 3, relative position of these shapes is also optimized by 2D/3D registration proposed in Section 2.3, and Average error 3 becomes
smaller than Average error 2.

parameters in ﬂuoroscopic images using the proposed techniques.
Examples of the ﬂuoroscopic images are shown in Fig. 11.
Fig. 12 shows one example of the pair (No. 4 and 24) of ﬂuoroscopic images, which were captured from directions crossing atright angles. The average errors of the estimated femoral shapeare shown in Fig. 13 and Table 2 for various numbers of shapeparameters used for the estimation. The estimation process andthe estimated 3D shape when the number of the estimated shapeparameters is 10 are shown in Figs. 14 and 15. The calculation timeby a 3.2-GHz Pentium IV, is approximately one minute, including
the contour detection by the Level Set Method and the shapeparameter estimation.
Finally, the average errors for various pairs of ﬂuoroscopic
images are shown in Fig. 16 when the number of the estimatedshape parameters is 10. As a result of a series of experiments using
the phantom femur, we concluded that the 3D shape can be estimated with an average error of less than 1.2 mm if we choose
the proper images captured from directions that intersect approximately at right angles.
3.3. In vivo experiments
Next, we conducted in vivo experiments for hip prosthesis patients. Fluoroscopic images of four patients were taken in clinicalpractice using the ﬂuoroscopic imaging apparatus (Siemens, Sire-

3.2. Experiments using the phantom femur
We conducted experiments using a dry femur bone and ﬂuoroscopic images. In these experiments, a special ﬂuoroscopic imaging
apparatus (Siemens, Siremobil ISO-C) was used for the ﬂuoroscopic
photography from various directions around the phantom femur.
First, we captured images of calibration markers with nine glass
bubbles (left of Fig. 10) at 50 positions around the markers from 0
to 190° using the ﬂuoroscopic apparatus. The 3D positions of themarkers were also measured precisely by the CT scanner. Next,
the intrinsic and extrinsic parameters of the ﬂuoroscopic apparatus
were calibrated by Tsai’s method [35].
After calibration, we replaced the markers with the dry femur
bone and captured 50 images at the same positions. In addition,
the precise 3D shape of the dry bone was measured by the CT scanner. Next, we chose two ﬂuoroscopic images from 50 images, asmentioned below, and estimated the position and the optimum

No.1
0 deg.

No.25
95 deg.
(approx.)Femoral dry bone
No.50
190 deg.

Markers

Calibration markers

Fig. 11. Measured ﬂuoroscopic images.

No.4
Fig. 10. Fluoroscopic images of the calibration markers and the phantom femur.

Femur

No.24

Fig. 12. Examples of ﬂuoroscopic images (Nos. 4 and 24).

208

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211

2.5

[mm]

Average error [mm]

14

24

34

44

1.246

0.992

1.128

1.377

0.980

1.050

1.223

1.143

1.091

4

3D CT image(Average error 2)

2

Image No.

1.5

4 (∼ 15 deg.)1

two radiographs(Average error 3, Proposed method)

0.5

AS : Average shape
0

AS

1

2

3

4

5

6

7

8

9

10

14 (∼ 53 deg.)-

24 (∼ 91 deg.)--

34 (∼ 129 deg.)

-

-

-

44 (∼ 167 deg.)

-

-

-

1.504
-

Number of parameters
Fig. 16. Estimation errors for various pairs of ﬂuoroscopic images.
Fig. 13. Comparison of estimation errors using CT image and two radiographs.

mobil ISO-C), and shapes measured by the CT scanner and estimated by parameter estimation were compared.
Table 2Average of estimation errors of the phantom femur (mm)Number of parameters

Average

STD.

Max.

0
1
2
3
4
5
6
7
8
9
10

1.544
1.568
1.344
1.103
1.068
1.020
0.986
1.009
0.929
0.931
0.889

1.009
1.021
1.111
0.868
0.869
0.814
0.817
0.785
0.743
0.738
0.709

5.381
5.486
6.653
4.820
4.910
4.776
4.770
4.174
3.813
3.709
3.758

Fig. 14. Process of shape and position estimation of the femur.

3.3.1. Calibration procedure
First, we measured the internal and external parameters of the
ﬂuoroscopic apparatus. The internal parameters consists of ﬁveparameters, that is, focus length, 1st order radial lens distortion,center pixel x and y, and scale factor. On the other hand, the external parameters contain six parameters, that is, the position and orientation of a radiation source in 3D space. After the calibration
process, we can get correct ﬂuoroscopic images of known scale
without lens distortion, and the relative positions of two viewpoints. These parameters are all requisite for the proposed shapeestimation procedure. A non-coplanar marker and a coplanar marker, shown in Fig. 17, are used for the calibration of external andinternal parameters, respectively. These markers are constructedof acrylonitrile butadiene styrene (ABS), which indicated the highest transmission of X-rays in preliminary experimentation. The
non-coplanar marker contains nine small stainless steel spheres
and the coplanar marker contains 16 small stainless steel disks.
The calibration procedure using these markers is as follows:
Step 1. The non-coplanar marker is captured by X-ray CT and therelative positions of the stainless steel spheres are measured.
Step 2. A ﬂuoroscopic image of the non-coplanar marker is captured, as shown in Fig. 18, and the internal parameters of
the ﬂuoroscopic imaging apparatus are estimated by Tsai’s
method.
Step 3. The coplanar marker is placed under the patient’s hip andtwo ﬂuoroscopic images of the femur are captured fromtwo directions, as shown in Figs. 19 and 20.

ABS
90 mm
50 mm

150 mm

ABS

Stainless steel ball
Stainless steel disk

Non-coplanar marker for
external parameter calibration
Fig. 15. Actual femoral shape and estimated femoral model with and without shapeparameter estimation.

Coplanar marker forinternal parameter calibration

Fig. 17. Non-coplanar and coplanar markers.

209

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211

Fig. 21. Extracted contours of the femur in ﬂuoroscopic images.

Fig. 18. Example ﬂuoroscopic image of non-coplanar marker.

X-ray apparatus

position and 10 shape parameters from the silhouette of the femur.
The precise 3D shapes of patients’ femurs were measured precisely
by CT scanner.
The average errors and standard deviations of the estimatedfemoral shapes are shown in Fig. 22. In these ﬁgures, ‘‘A” in thehorizontal axis shows the case in which only the position is estimated without parameter estimation.
The experimental results show that the average error between
the estimated shape and the actual shape is approximately from
0.8 to 1.1 mm for the in vivo experiments. One example of average,
actual, and estimated shapes for Case 4 is shown in Fig. 23, and the
distribution of average error is shown in Fig. 24. In Fig. 24, dark regions indicate less error, and the brightness of each point is proportional to its average error. From this ﬁgure, we veriﬁed that theerrors in the femoral head and lesser trochanter are reduced.

Coplanar marker

1.6

case1
1.4

Average error [mm]case2Fig. 19. X-ray photography of coplanar marker.

1.2
1.0

case3

0.8

case4

0.6
0.4
0.2

AS : Average shape
0
AS

1

2

3

4

5

6

7

8

9

10

11

12

11

12

Number of parameters
1.2

case1

case2

1.0
Fig. 20. Two ﬂuoroscopic femoral images of hip prosthesis patients.

S.D. [mm]0.8Step 4. The external parameters of two ﬂuoroscopic images are
estimated using the internal parameters obtained in Step
2 and the projection image of stainless steel disks by Tsai’s
method. After this calibration process, the relative positions of two viewpoints of ﬂuoroscopic images are precisely
determined.

case3

0.6

case4

0.4
0.2

AS : Average shape
3.3.2. Results
First, we manually extracted the contour of the femur (Fig. 21)
in the ﬂuoroscopic images (Fig. 20). Image size and resolution are
512  512 pixels and 0.454 mm pixel. We then estimated the

0
AS

1

2

3

4

5

6

7

8

9

10

Number of parameters
Fig. 22. Average error and standard deviation for the number of estimated shapeparameters for the femurs of four patients.

210

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211

was veriﬁed that the 3D shape can be estimated using two 2D ﬂuoroscopic images taken from different view points with the sameaccuracy as in the case of 3D shapes being compared directly.
In the future, optimum conditions, such as the optimum number and directions of the ﬂuoroscopic images, will be investigated,and clinical experiments in ﬂuoroscopic image diagnosis will be
performed.
Appendix A
A.1. Fast extraction of the projected contour line of the 3D parametricfemoral modelFig. 23. Actual and estimated femoral model with and without shape parameter
estimation for Case 4.

In Step (4) in Section 2.3.1, the contour detection and identiﬁcation triangular patches on the 3D model corresponding to points on
the contour line are computationally expensive and time consuming. In our implementation, we utilize the high-speed renderingfunction of the OpenGL hardware accelerator, and thus these procedures are executed quite rapidly.
The detailed algorithm is as follows. Initially, we assigned different colors to all of the triangular patches in the 3D model anddraw the projected image of the 3D model on the image bufferusing the OpenGL hardware accelerator. The contour points of
the 3D model are detected by raster scanning of the image buffer.
By reading the colors of the detected contour points, we can identify the corresponding triangular patches on the 3D geometricmodel.
References

Fig. 24. Error distribution indicated by brightness (higher brightness indicateslarger error.

4. Conclusions
We proposed a method by which to estimate the 3D shape of
the in vivo femur from only two ﬂuoroscopic images using a parametric femoral model. Although the precise 3D shape of the femur
is usually measured using a CT scanner or an MRI system, the proposed method enables a precise 3D shape to be estimated using
only two ﬂuoroscopic images taken by an inexpensive ﬂuoroscopic
inspection apparatus. Thus, the cost of the inspection system canbe dramatically reduced and the 3D image-based medical diagnosis becomes available even in small clinics.
In vivo experiments revealed the average error between theestimated shape and the actual shape to be 0.8–1.1 mm, and it[1] I. Stamos, P.K. Allen, Integration of range and image sensing for photorealistic
3D modeling, in: Proceedings of the 2000 IEEE International Conference on
Robotics and Automation, 2000, pp. 1435–1440.
[2] I. Stamos, P.K. Allen, Automatic registration of 2D with 3D imagery in urban
environments, in: Proceedings of the International Conference on Computer
Vision, 2001, pp. 731–737.
[3] L. Liu, I. Stamos, Automatic 3D to 2D registration for the photorealistic
rendering of urban scenes, in: IEEE International Conference on Robotics andAutomation, 2005.
[4] R. Kurazume, K. Noshino, Z. Zhang, K. Ikeuchi, Simultaneous 2D images and 3D
geometric model registration for texture mapping utilizing reﬂectanceattribute, in: Proceedings of Fifth Asian Conference on Computer Vision(ACCV), 2002, pp. 99–106.
[5] M.D. Elstrom, P.W. Smith, Stereo-based registration of multi-sensor imagery
for enhanced visualization of remote environments, in: Proceedings of the1999 IEEE International Conference on Robotics and Automation, 1999, 1948–
1953.
[6] K. Umeda, G. Godin, M. Rioux, Registration of range and color images using
gradient constraints and range intensity images, in: Proceedings ofSeventeenth International Conference on Pattern Recognition, 2004, 12–15.
[7] Y. Iwashita, R. Kurazume, K. Hara, T. Hasegawa, Fast alignment of 3D
geometrical models and 2D color images using 2D distance maps, in:Proceedings The Fifth International Conference on 3D Digital Imaging andModeling (3DIM), 2005, pp. 164–171.
[8] Q. Delamarre, O. Faugeras, 3D articulated models and multi-view tracking withsilhouettes, in: Proceedings of the International Conference on Computer
Vision, vol. 2, 1999, pp. 716–721.
[9] K. Matsushita, T. Kaneko, Efﬁcient and handy texture mapping on 3D surfaces,
In Comput. Graphics Forum 18 (1999) 349–358.
[10] P.J. Neugebauer, K. Klein, Texturing 3D models of real world objects frommultiple unregistered photographic views, In Comput. Graphics Forum. 18(1999) 245–256.
[11] G. Penny, J. Weese, J. Little, P. Desmedt, D. Hill, D. Hawkes, A comparison ofsimilarity measures for use in 2-D/3-D medical image registration, IEEE Trans.
Med. Imaging 17 (4) (1998) 586–595.
[12] L. Zollei, E. Grimson, A. Norbash, W. Wells, 2D–3D rigid registration of X-ray
ﬂuoroscopy and CT images using mutual information and sparsely sampled
histogram estimators, in: Proceedings of Computer Vision and Pattern
Recognition, 2001, pp. 696–703.
[13] C.V. Stewart, C.L. Tsai, A. Perera, A view-based approach to registration:
theory and application to vascular image regisrtaion, in: International
Conference on Information Processing in Medical Imaging (IPMI), 2003, pp.
475–486.
[14] A. Guimond, A. Roche, M. Ayache, J. Meunier, Three-dimensional multimodal
brain warping using the demons algorithm and adaptive intensity corrections,
IEEE Trans. Med. Imaging 20 (1) (2001) 58–69.

R. Kurazume et al. / Computer Vision and Image Understanding 113 (2009) 202–211[15] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality
image registration by maximization of mutual infomation, IEEE Trans. Med.
Imaging 16 (2) (1997) 187–198.
[16] C.V. Stewart, C.L. Tsai, A. Perera, Rigid and afﬁne registration of smoothsurfaces using differential properties, in: Proceedings of Third European
Conference on Computer Vision (ECCV’94), 1994, pp. 397–406.
[17] C.R. Meyer, J.L. Boes, B. Kim, P.H. Bland, K.R. Zasadny, P.V. Kison, K. Korak, K.A.
Fery, R.L. Wahl, Demonstration of accuracy and clinical versatility of mutual
information for automatic multimodality image fusion using afﬁne and thinplate spline warped geometric deformations, Med. Image Anal. 1 (3) (1997)
195–206.
[18] A. Gueziec, X. Pennec, N. Ayache, Medical image registration using geometrichashing, IEEE Comput. Sci. Eng. Spec. Issue Geometric Hashing 4 (4) (1997) 29–
41.
[19] P.R. Andresen, M. Nielsen, Non-rigid registration by geometry constrained
diffusion, Med. Image Comput. Comput. Assist. Interv. (MICCAI’ 99) (1999)
533–543.
[20] T. Masuda, Y. Hirota, K. Ikeuchi, K. Nishino, Simultaneous determination ofregistration and deformation parameters among 3D range images, in: Fifth
International Conference on 3-D Digital Imaging and Modeling, 2005, 369–376.
[21] C.S.K. Chan, D.C. Barratt, P.J. Edwards, G.P. Penney, M. Slomczykowski, T.J.
Charter, D.J. Hawkes, Cadaver validation of the use of ultrasound for 3D model
instantiation of bony anatomy in image guided orthopaedic surgery, in:
Lecture Notes in Computer Science, 3217 (Proceedings Seventh International
Conference on Medical Image Computing and Computer Assisted Intervention,
Part II (MICCAI 2004), St.-Malo, France), 2004, pp. 397–404.
[22] D. Terzopoulos, D. Metaxas, Dynamic 3D models with local and global
deformations: deformable superquadrics, IEEE Trans. Pattern Anal. Mach.
Intell. 13 (7) (1991) 703–714.
[23] K. Nakamura, R. Kurazume, T. Okada, Y. Sato, N. Sugano, T. Hasegawa, 3D
reconstruction of a femoral shape using a parametric model and two 2D
radiographs, in: Proceedings of Meeting on Image Recognition andUnderstanding, 2006, pp. 78–83.
[24] G. Zheng, M. Ballester, M. Styner, L. Nolte, Reconstruction of patient-speciﬁc
3D bone surface from 2D calibrated ﬂuoroscopic images and point distribution

211

model, in: Proceedings of Medical Image Computing and Computer-Assisted
Intervention (MICCAI’06), 2006, pp. 25–32.
[25] Keita Yokota, Toshiyuki Okada, Masahiko Nakamoto, Yoshinobu Sato,
Masatoshi Hori, J. Masumoto, Hironobu Nakamura, Shinichi Tamura,
Construction of conditional statistical atlases of the liver based on spatial
normalization using surrounding structures, Int. J. Comput. Assist. Radiol. Surg.
(2006) 39–40.
[26] Toshiyuki Okada, Ryuji Shimada, Yoshinobu Sato, Masatoshi Hori, Keita
Yokota, Masahiko Nakamoto, Yen-Wei Chen, Hironobu Nakamura, Shinichi
Tamura, Automated segmentation of the liver from 3D CT images using
probabilistic atlas and multi-level statistical shape model, in: Proceedings of
the Tenth International Conference on Medical Image Computing andComputer Assisted Intervention (Proc. MICCAI 2007), 2007, pp. 86–93.
[27] T.F. Cootes TF, C.J. Cooper, C.J. Taylor, J. Graham, Active shape models—theirtraining and application, Comput. Vis. Image Understanding 61 (1) (1995) 38–
59.
[28] T. Okada, M. Nakamoto, Y. Sato, N. Sugano, H. Yoshikawa, S. Tamura, T. Asaka,
Y.-W. Chen, Effects of surface correspondence methods in statistical shapemodelling of the proximal femur on approximation accuracy, in: Proceedingsof the Twenth International Congress and Exhibition Computer Assisted
Radiology and Surgery CARS, 2006, p. 016.
[29] Haili Chui, Anand Rangarajan, A new point matching algorithm for non-rigid
registration, Comput. Vis. Image Understanding 89 (2-3) (2003) 114–141.
[30] M. Greenspan, M. Yurick, Approximate k-d tree search for efﬁcient ICP, in:Proceedings of 3-D Digital Imaging and Modeling 3DIM, 2003, pp. 442–448.
[31] J. Sethian, Level Set Methods and Fast Marching Methods, second ed.,Cambridge University Press, UK, 1999.
[32] J. Sethian, A fast marching level set method for monotonically advancing fronts,
in: Proceedings of the National Academy of Science 93, 1996, pp. 1591–1595.
[33] G. Borgefors, Distance transformations in arbitrary dimensions, Comput. Vis.
Graphics Image Process. 27 (3) (1984) 321–345.
[34] S.W. Shih, Y.T. Wu, Fast Euclidean distance transformation in two scans using a
3 3 neighborhood, Comput. Vis. Image Understanding 93 (2) (2004) 195–205.
[35] R.Y. Tsai, An Efﬁcient and Accurate Camera Calibration Technique for 3D
Machine Vision, in: Proceedings of IEEE Conference on Computer Vision and
Pattern Recognition, 1986, pp. 364–374.