Automatic Extraction of Femur Contours from ... - Semantic Scholar

46

JOURNAL OF MULTIMEDIA, VOL. 2, NO. 5, SEPTEMBER 2007

Automatic Extraction of Femur Contours from Calibrated X-Ray Images using Statistical Information Xiao Dong, Miguel A. Gonzalez Ballester, Guoyan Zheng MEM Research Center, Bern University, Bern, Switzerland Email: {guoyan.zheng, xiao.dong}@MEMcenter.unibe.ch

Abstract— Automatic identification and extraction of bone contours from x-ray images is an essential first step task for further medical image analysis. In this paper we propose a 3D statistical model based framework for the proximal femur contour extraction from calibrated x-ray images. The automatic initialization to align the 3D model with the x-ray images is solved by an Estimation of Bayesian Network Algorithm to fit a simplified multiple component geometrical model of the proximal femur to the x-ray data. Landmarks can be extracted from the geometrical model for the initialization of the 3D statistical model. The contour extraction is then accomplished by a joint registration and segmentation procedure. We iteratively updates the extracted bone contours and an instanced 3D model to fit the x-ray images. Taking the projected silhouettes of the instanced 3D model on the registered x-ray images as templates, bone contours can be extracted by a graphical model based Bayesian inference. The 3D model can then be updated by a non-rigid 2D/3D registration between the 3D statistical model and the extracted bone contours. Preliminary experiments on clinical data sets verified its validity. Index Terms— contour extraction, statistical model, Bayesian network, 2D/3D registration, segmentation, calibrated x-ray image

I. I NTRODUCTION X-ray images are still playing a crucial role in diagnosis and surgery. Accurate extraction of bone contours from x-ray images is an essential component of the computer analysis of medical images for diagnosis [1] [2] [3], planning [4] [5] [6] or 3D reconstruction of anatomic structures [7] [8] [9]. X-ray images can vary a lot in the image brightness and contrast as well as in the imaged region of anatomy. The homogeneity assumption on the intensity distribution, on which most of the segmentation algorithms depend, are not satisfied for large and complex bones in x-ray images. Occlusion such as the overlap between the femoral head and the pelvis bone also happens so that local features of bone contours can not be reliably identified. Therefore conventional segmentation techniques [1] [5] [6], which mainly depend on local image features such as the edge information, can not provide a satisfactory solution. Model based segmentation A part of the work has been presented at and published in the proceedings of IEEE Workshop on Applications of Computer Vision (WACV) 2007, February 21 - 22, 2007, Austin, Texas, USA

© 2007 ACADEMY PUBLISHER

techniques are usually implemented to obtain robust and accurate results [3] [7] [10] [11] [12]. In [3] [11] [12] [13] [14], 2D statistical models are constructed from a training image set under the assumption that the images are taken from a certain view direction. 2D statistical models can encode both the shape and texture information learnt from training data set, which is helpful to improve the robustness and accuracy of the segmentation of noisy images. But the 2D statistical model asks for a proper initialization due to its limited convergence region. Fully automatic initialization can be accomplished by the generalized Hough transformation [11], neural nets [12] or evolutionary algorithms [13] [14] [15]. But both the initialization and segmentation performance rely on whether the view direction assumption can be fulfilled. 3D statistical models [7] [8] [9] have been used for 2D segmentation and 3D reconstruction from calibrated 2D x-ray images. The advantage of their application in segmenting an image taken from an arbitrary view direction is apparent.3D statistical model also need an initialization, which is usually achieved by user interactions in current systems [7] [8] [16]. Due to the dense mesh of the 3D statistical model, fully automated initialization based on evolutionary algorithms is computational expensive [15]. In this paper we propose a 3D statistical model based fully automatic segmentation framework for the proximal femur bone contour extraction from calibrated x-ray images. The automatic initialization is accomplished by first fitting a simplified multiple component geometrical model of the proximal femur to the x-ray images and then aligning the 3D statistical model to fit the geometrical model. The statistical 3D model based contour extraction is then solved by a joint segmentation/non-rigid registration procedure, which iteratively updates the extracted contours and the instanced 3D model to fit the x-ray images. In our work we use the projected silhouette of the 3D statistical model on each x-ray image as a 2D template and accomplish the template based contour extraction by a Bayesian inference on a graphical model. The update of the 3D model is achieved by a 2D/3D registration to fit an instanced 3D model with the extracted contours. The remainder of this paper is organized as follows: In Section 2, related work are categorized and discussed. Then, in Section 3, we give an introduction of the calibrated x-ray images and the Principle Component


Analysis (PCA) based statistical model. In Section 4 the Estimation of Bayesian Network Algorithm (EBNA) for the automatic initialization of the statistical model is given and the 3D statistical model assisted bone contour extraction is described in Section 5. Results of experiments on clinical data set are discussed in Section 6. Finally we conclude our work in Section 7. II. R ELATED W ORK The problems we are dealing with are essentially MAP estimations, for which Bayesian inference on graphical models works as a general framework. To adopt the framework for a specific application, the task is then to select the proper graphical model to represent the stochastic system and to design an efficient Bayesian inference algorithm to accomplish the MAP estimation. Various algorithms, which can all be interpreted in the language of Bayesian inference on graphical models including particle filtering [17], particle smoothing [18], belief propagation (BP), loopy belief propagation (LBP) and generalized belief propagation (GBP) [19], were invented or adopted in different research fields. In statistical physics, replica Monte Carlo [20] was used on a k-dimensional grid lattice to simulate the phase transition or to find the ground state of the spin glass. In computer science field, survey propagation algorithm [21] was proposed to solve the K-satisfiability problem and belief propagation on a junction-tree was used for graph matching [22]. In signal processing field, a Markov chain was established for the state-space model of the sequentially received signal, on which particle filtering and particle smoothing were implemented to calculate different marginal posterior distributions to solve the model selection problem [17] or to estimate the direction of arrival(DOA) in array signal processing [23]. In telecommunication signal processing, belief propagation and generalized belief propagation on factor graphs were discussed in the decoding of error-correcting codes such as Turbo-code [24] and low density parity check code(LDPC) [25]. In control theory, assumed-density filtering (ADF) and expection propagation (EP) [26] were proposed to approximate posteriors in Bayesian networks. In artificial intelligence or machine learning, Estimation of Distribution Algorithm(EDA) [27] and Estimation of Bayesian Network Algorithm(EBNA) were implemented on the unsupervised learning of Bayeisan network. In computer vision fields, CONDENSATION algorithm [28], essentially a particle filter, worked on a Markov chain graph of the state space model of the video frames to track object contours. In [29] [30] [31] [32] a multiple component object (human body, hand) with its temporal behavior was model by a dynamic Bayesian network (DBN). Nonparametric belief propagation (NBP) [31] and sequential mean-field Monte Carlo (MFMC) [30] [33] were implemented to accomplish the posteriori probability inference in continuous valued state space for object tracking. Similar to non-parametric belief propagation, PAMPAS [34] is an alternative for continuous valued probability inference method with combines © 2007 ACADEMY PUBLISHER

47

belief propagation with the idea from particle filtering. A data-driven Markov Chain Monte Carlo(MCMC) was used in [35] to estimate human upper body pose in static images. In [3] the lumbar vertebra were segmented by linked Active Appearance Model (ASM). This work also involves information passing among components and can be regarded as a smoothing algorithm for a state-space model. In image processing field, Bayesian network is also exploited for finding deformable shapes [36] [37], where both the local feature and global shape information can be encoded in a graphical model and belief propagation is carried out to find the solution. III. C ALIBRATED X - RAY IMAGES AND THE STATISTICAL MODEL

A. Calibrated x-ray images The calibrated x-ray image is an extension of the traditional x-ray image so that not only the x-ray image but also the projection parameters of the x-ray shot are recorded [16]. Therefore the calibrated x-ray images are co-registered with respect to a common reference coordinate system. In our experiment, calibrated x-ray images from C-arm are used. Due to the limited imaging volume of C-arm, we ask for four images of the proximal femur from different view directions, of which two images focus on the femoral head and the other two focus on the femoral shaft. The calibrated x-ray image set is noted by I = {Ii , Pi }i=1,2,3,4 , where Ii is the ith x-ray image and Pi is its correspondent projection parameters. B. Statistical model of the proximal femur A Principle Component Analysis (PCA) based 3D statistical sufrace model M with 4098 vertices of the proximal femur is constructed from a training data set containing the CT data of 13 bones [16] as shown in Fig. 1. An instance generated from the statistical model with parameter set Q = {α, β0 , β1 , . . . , β11 } can be described as 11 X 1 (1) M : S(Q) = α(S0 + βi λi2 Pi ) i=0

where S0 = [S00 , . . . , S04097 ]T is the mean shape of the model, α is the scaling factor, λi and Pi are the ith eigenvalue and the the correspondent eigenvector of the correlation matrix of the training data. IV. AUTOMATED INITIALIZATION OF THE STATISTICAL MODEL

To find the initial rigid transformation T0 and parameter set Q0 to align the model instance S(Q0 ) with the observed x-ray images, a multiple component geometrical model is constructed for the approximal femur. A Bayesian network is established to encode the constraints among the components and an Estimation of Bayesian Network Algorithm (EBNA) to align the geometrical model with the x-ray images. Then T0 and Q0 can be calculated from the geometrical model accordingly.

48


anatomical structure of the proximal femur. In our approach they are all set to be uniform distributions in their allowed parameter spaces based on anatomy knowledge. C. Geometrical model fitting by EBNA Fitting the geometrical model with the x-ray images can be easily interpreted as a MAP estimation (X∗geo ) = max P rob(Xgeo |I) Xgeo

Figure 1. PCA based 3D statistical surface model of the proximal femur

(2)

where taking the priori distribution of Xgeo as a uniform distribution, we have P rob(Xgeo |I) ∝ P rob(I|Xgeo ) and the latter is called the observation model. 1) Observation model: The observation model is constructed by a similarity measurement between the x-ray images and the silhouettes of the projected geometrical model on the correspondent image planes. It can be computed on sampled positions along the projected silhouettes in a similar way as the similarity measurement of Adaptive Appearance Model (AAM) as follows P rob(I|Xgeo ) =

4 Y

pE (Ii |Xgeo )pI (Ii |Xgeo )

(3)

i=1

(a) Multiple component geometrical model for the proximal femur

where pE (Ii |Xgeo ) and pI (Ii |Xgeo ) are the edge matching likelihood and the intensity matching likelihood between the ith observed x-ray image and the projected silhouette of the geometrical model on the ith x-ray image respectively. pE (Ii |Xgeo ) = exp(−λE

X

(dE (ekp , eIi ))2 )

(4)

k

(b) Graphical model for the geometrical model fitting Figure 2. Models for automatic initialization

where dE (ekp , eIi ) is the shortest distance between a sampled position ekp on the projected silhouette of the projected geometrical model and the edge point set of the ith x-ray image eIi . λE is a control parameter. And similarly X (5) pI (Ii |Xgeo ) = exp(−λI (dI (ekp , Ii ))2 ) k

A. Multiple component geometrical model of the proximal femur The proximal femur is modeled by a geometrical model consisting of 3 components: head, neck and shaft, which are described by a sphere, a trunked cone and a cylinder with parameter set Xgeo = {XH , XN , XS } respectively as shown in Fig. 2(a). The configurations of the three components are constrained by the anatomical structure of the proximal femur. B. Bayesian network for the geometrical proximal femur model The constraints among components are encoded in the potentials among the nodes in a graphical model [35] [31] [30] as shown in Fig. 2(b). The potentials are set to constrain the distances and angles between components so that the geometrical model can represent a meaningful © 2007 ACADEMY PUBLISHER

dI (ekp , Ii )

is a metric by which the similarity where between the intensity distribution of the simulated xray image(by projecting the current geometrical model on the image plane using the projection parameters of each x-ray shot and assuming the intensity distributions of the interior/exterior of the silhouette are dark/bright with fixed intensity values) at the neighborhood of the sampled position ekp and the intensity distribution at the neighborhood of the same position on the x-ray image Ii . Belief propagation based approaches such as nonparametric belief propagation and sequential mean field Monte Carlo [30] [31] ask for a proper initialization so that most of the particles are located near the global optimal solution. Therefore most of the messages carry meaningful and efficient information among nodes of the graphical model to help the algorithm converge. Without a proper initialization, the message passing will show a very low efficiency. Data driven MCMC [35] uses low level feature


49

detectors to identify component candidates of the object from the data to generate proposal distributions for particles to improve the efficiency. But its performance highly depends on the correctness of the low level detectors. In our problem to identify the proximal femur from x-ray images, no proper initialization is available and the detection of proximal femur components based on low-level detectors is not reliable due to the low image quality. We use an Estimation of Bayesian Network Algorithm (EBNA) to accomplish the MAP estimation as given in Algorithm 1. 1. Initialization Generate the first generation of particle set with N particles {X0geo,i }i=0,...,N −1 from the proposal distributions q 0 (Xgeo ) = q 0 (XS )p0 (XN |XS )p0 (XH |XN ) q 0 (XS ) = π(XS ) p0 (XN |XS ) = π(XN )ψ(XS , XN ) p0 (XH |XN ) = π(XH )ψ(XN , XH )

(a) Fitting the geometrical model with x-ray images

2. Observation Given the current generation of particle set, calculate the weight of each particle as win ∝ P rob(I|Xn geo,i ). 3. Update Update the proposal distributions as q n+1 (Xgeo ) = q n+1 (XS )pn+1 (XN |XS )pn+1 (XH |XN ) q n+1 (XS ) = NPDE(win , Xn S,i ) pn+1 (XN |XS ) = NPDE(win , Xn N,i )ψ(XS , XN ) pn+1 (XH |XN ) = NPDE(win , Xn H,i )ψ(XN , XH ) where NPDE(win , Xn S/N/H,i ) is a nonparametric density estimation [38]. Generate the next generation of particle set from the updated proposal distributions. 4. Go to 2 until the particle set converges. 5. After convergence, the configuration of the particle with the highest weight is selected as our final result.

Algorithm 1. EBNA for geometrical model fitting In fact, our EBNA can also be interpreted as a particle filter if we regard each iteration of the EBNA as a time slot of a dynamic Bayesian network, for which the temporal evolution of the state space is static and the observation of each slot remain the same as our observed x-ray images. The proposal method of the particle filter can then be regarded as a partitioned sampling as in [39]. The reason that the partition sequence of the state space as XS : XN : XH is that it’s more reliable to identify the femoral shaft than the femoral head due to the occlusion at hip joints. D. Statistical model initialization From the mean shape of the 3D statistical model S0 , the model vertices can be classified into three regions, femoral head, neck and shaft. The femoral head center and radius, axes of femoral neck and shaft can be determined in the model coordinate space by a 3D sphere fitting to the femoral head region and cylinder fittings to the femoral neck and shaft regions. From the geometrical model, landmarks including the femoral head center, femoral head radius, femoral neck axis and femoral shaft axis, can be easily extracted. On both the geometrical model and the mean shape of the statistical model, a reference © 2007 ACADEMY PUBLISHER

(b) Statistical model initialization by fitting it with the geometrical model Figure 3. Automatic 3D statistical model initialization

coordinate system can be established taking the femoral head center as the origin and the plane defined by the femoral neck and shaft axes as x-y plane. The initial rigid transformation T0 can then be computed to align the two coordinate systems. The scale factor α0 can be S G S G +DHC estimated by the ratio (RH SA )/(RH +DHC SA ). S G RH and RH are the femoral head radii of the mean shape of the statistical model and the fitted geometrical G S model. DHC SA and DHC SA are the distances between the femoral head center and the femoral shaft axis. The estimated α0 T0 can then roughly fit the statistical model(the scaled mean shape) to the geometrical model as shown in Fig. 3(b). V. 3D STATISTICAL MODEL BASED CONTOUR EXTRACTION

After the statistical model initialization, the contour extraction is accomplished by a joint registration and segmentation as summarized in Algorithm 2 and Figure 4. A. Simulated x-ray and silhouette extraction Given the current instanced 3D model M : S(Qn ) and the transformation Tn , a simulated x-ray can be generated by projecting each triangle patch on the mesh surface of the transformed 3D model on the image plane of each x-ray image Ii using the projection parameter Pi of that x-ray image as shown in Fig. 5(a). The silhouette of the simulated x-ray image can be easily obtained using a Canny edge detector as shown in Figure 5(b).

50


(a) Projected 3D model superimposed with x-ray images

(b) Extracted silhouette of the projected 3D model superimposed with the observed x-ray images, due to the limited femoral shaft length of the statistical model, the silhouettes of the femoral shaft are trunked to get rid of the undesired edge segment from the open end of the femoral shaft Figure 5. Simulated x-ray and silhouette extraction

Figure 4. Schematic plot of the 3D statistical model based contour extraction

1. Simulated x-ray and silhouette extraction Given the current instanced statistical model M : S(Qn ) and the transformation Tn to align the 3D model to the observed x-ray images, project the aligned statistical model on each of the K(K ≥ 2) x-ray image planes using the projection geometry of each x-ray image. From the simulated x-ray images the silhouettes {Ck,n model }k=0,...,K−1 are extracted [8]. 2. 2D template based segmentation On each x-ray image, taking the correspondent silhouette of the projected statistical model Ck,n model as a template, a Bayesian network based shape matching is implemented to search for the bone contour Ck,n image . 3. Nonrigid 2D/3D registration A 2D/3D nonrigid registration is carried out to fit the extracted bone contours Ck,n image k=0,...,K−1 and the statistical model M, which results in an updated instanced model M : S(Qn+1 ) and a rigid transformation Tn+1 . 4. Go to 1, until the procedure converges.

Algorithm 2. Statistical model based contour extraction B. 2D template based segmentation using belief propagation Usually active contour or statistical model (ASM or AAM) [1] [10] are used for model based 2D segmentation. Ideally the model should keep the global shape information as well as involve the local feather information such as the edge and intensity distribution. Active contour emphasizes on searching for local features along © 2007 ACADEMY PUBLISHER

the contour and therefore usually lacks of the ability to hold the global structure. ASM/AAM consider both the global shape and the local features. But its shortcoming is that its searching strategy usually independently estimate each point to a new target position by a local search along the normal direction of the contour. Therefore the inaccuracies in this estimation can not be accommodated for during the local search by the global shape information and can only be regularized by a projection to the shape space during the nonrigid registration step. In [36] [37] the shape matching problem is formalized as an Bayesian inference on a graphical model and solved by loopy belief propagation and Bethe free energy approximation respectively. In this approach, a graphical model is established so that the correspondence assignment for each point involves both the global shape and local feature information. From the silhouette of the projected 3D statistical model, we sample M points(nodes) tracing along the contour as the shape prior. Each point is described by a parameter set qi = {xi , gi , fi }, i = 1, . . . , M , where xi = (xi , yi ) is the position of ith point in the image coordinate system, gi = (gxi , gyi ) is the gradient vector of the x-ray image,fi = 1 if the current node belongs to the femur head projection silhouette and fi = 0 otherwise. The configuration of our model can then be written as Qmodel = {qi }i=1,...,M , where gi is set as the tangent vector of the template curve at position xi . The configuration of a candidate contour can be written 0 as Qcand = {qi }i=1,...,M . We then establish a partially connected graph with M vertices as: G(V, E), V =


51

{vi , }i=1,...,M , E = {eij }i,j=1,...,M ,where  1, |i − j| ≤ NShaf t , fi = 0, fj = 0, i 6= j    1, f = 0, f = 1 i j eij = (6)  1, |i − j| ≤ NHead , fi = 1, fj = 1, i 6= j    0, other NHead and NShaf t determine the number of connected neighbors for the head nodes and the non-head nodes respectively. Larger NHead and NShaf t will keep the rigidity of the shape but will fail to track deviation from the template. The reason that all the head nodes are connected with the non-head nodes is that we need the non-head nodes(which are supposed to be relatively easier to be located than the head nodes) to guide the localization of the head nodes. On each vertex, a local observation ψi (qi ) is defined and for each eij = 1 a pairwise potential is defined as ψij (qi , qj ). The correspondent factor graph [19] is shown in Fig. 6(b). Given the template Qmodel = {qi }i=1,...,M , the joint probability distribution of the factor graph with an can0 didate configuration Qcand = {qi }i=1,...,M is then given by p(Qcand ) =

Y 0 0 0 1 Y ψi (qi ) ψij (qi , qj ) Z i e =1

(a) Graphical model built on the contour template

(b) Correspondent factor graph of the graphical model

(7)

ij

0

0

where ψ(qi ) = dot(gi , gi ),which means to penalize candidates with weak gradient amplitude and inconsis0 0 tent gradient direction with the model. ψij (qi , qj ) = −(µ

0 0 0 0 kkxi −xj k−kxi −xj kk (xi −xj )·(xi −xj ) +ν ) 0 0 kx i −xj k kx −x kkxi −xj k j i

which is set so e that the global shape of the model will be kept by penalizing the deviation of the angle and distance between vertices from our model. Under these definitions, a bone contour that keeps the global shape of our model and at the same time locates itself to the strong edge positions can be obtained by a MAP estimation as ∗ = Cimage

=

max

0

P (Qcand |I)

(8)

Qcand ={qi }

max

0 Qcand ={qi }

Y

0

ψi (qi )

i

Y

0

0

ψij (qi , qj )

(9)

eij=1

In our approach, the candidate positions for each node of the bone contour are sampled along the normal direction of the model and the standard loopy belief propagation [36] is used to approximate the MAP estimation and an example of the contour extraction using Bayesian inference is shown in Fig. 6(c).

C. 2D/3D nonrigid registration Given the extracted bone contours, the 3D statistical model can be updated by fitting it to the extracted bone k,n contours {Cimage } by a 2D/3D nonrigid registration procedure as described in Algorithm 3. An example of the nonrigid registration is shown in Fig. 7. © 2007 ACADEMY PUBLISHER

(c) Bayesian network based 2D segmentation, circles show the projected silhouettes and dots show the extracted contours Figure 6. Bayesian Inference on a graphical model for the contour extraction

1. Correspondence assignment For each point Pl on the extracted bone contour and the current projected silhouette of the current instanced statistical model M : S(Qn ) with its transformation Tn , a correspondence on the 2D image plane can be established by searching along the normal direction of the contour. Accordingly the correspondence between the back-projection line BP (Pl ) of Pl and a vertex vcorr(Pl ) on the model surface can be established as (BP (Pl ), vcorr(Pl ) ). Project vcorr(Pl ) on BP (Pl ) will generate a correspondent 3D point pair (vcorr(Pl ) , P roj(vcorr(Pl ) , BP (Pl ))). 2. Rigid transformation n+1 A rigid transformation Tupdate can be calculated from the paired point set {(vcorr(Pl ) , P roj(vcorr(Pl ) , BP (Pl ))),i=1,...,M } to align the current statistical model M : S(Qn ) to the extracted contours. The rigid transformation can then be n updated as Tn+1 = Tn+1 update T . 3. Constrained deformation of the 3D statistical model The residual error between correspondent point pairs can then be compensated by the constrained deformation of the statistical model [16] to find a updated 3D model M : S(Qn+1 ). 4. Go to 1, until the procedure converges.

Algorithm 3. 2D/3D non-rigid registration algorithm

52


TABLE I. S TATISTICAL RESULTS OF THE AUTOMATIC INITIALIZATION ALGORITHM , ALL RESULTS ARE RELATIVE TO THE MEAN VALUES OF THE

Figure 7. 2D/3D nonrigid registration to align the 3D model with the extracted bone contours

VI. E XPERIMENTAL R ESULTS We verified our approach on three set of clinical data, each data set includes four calibrated x-ray images of the proximal femur. The projection parameter of each x-ray shot is recorded by a projection matrix. It can be observed that in each data set the angle between the view directions of the two images focusing on the femoral head is quite small (less than 45 degrees), which will definitely increase the difficulty of identifying the 3D pose of the proximal femur. To test the robustness of the automatic initialization algorithm, we run the initialization algorithm for 10 trails on each data set using all the 4 images with particle number N = 200. In each trial the proximal femur is correctly identified and the statistical results are shown in Table I. The statistical results show a very robust performance of the initialization algorithm. We believe that the small variation among the results of different trials will have little influence on the following up contour extraction procedure. We then run our statistical model based contour extraction algorithm on the 2 images that focuse on the femoral head part in each data set with parameters set as M = 35, NHead = 4, NShaf t = 3, µ = 2, ν = 1 and the extracted proximal femur contours are shown in Fig. 8.

Parameter Head Center (mm) Head Radius (mm) Neck Radius (mm) Neck Length (mm) Neck Axis (degree) Shaft Radius(mm) Shaft Length(mm) Neck/Shaft Angle(degree)

10 TRIALS

Data Set 1 1.4±1.1 0.3±0.4 0.8±1.1 1.0±1.4 0.8±0.7 0.2±0.3 0.5±0.2 0.8±1.0

Data Set 2 0.1±0.1 0.6±0.2 0.6±0.9 1.3±1.8 2.3±1.0 0.1±0.2 0.9±0.5 2.0±2.5

Data Set 3 0.1±0.2 1.0±0.8 1.0±1.2 1.2±1.7 1.8±1.1 0.2±0.2 1.8±1.0 1.8±2.6

(a) Data set 1

(b) Data set 2

(c) Data set 3

VII. C ONCLUSIONS In this paper we proposed a 3D statistical model based fully automatic bone contour extraction framework from calibrated x-ray images. Our approach divides the MAP problem into two sub-problems, initialization and model based contour extraction. On each of them the graphical model based Bayesian inference plays a key role. We solve the automatic initialization by fitting a simplified multiple component geometrical 3D model to the observed x-ray images. The 3D model based initialization algorithm does not ask for strict view direction assumption compared with 2D model or 2D image feature based initialization. Since the fitting is accomplished by an evolutionary algorithm, it has a strong capability to overcome local optima and converge to the global optimum. The 3D statistical model based bone contour extraction is solved as a simultaneous 2D/3D registration and segmentation, where the template based segmentation is accomplished by a Bayesian inference procedure which in principle can outperform the active contour and © 2007 ACADEMY PUBLISHER

Figure 8. Results of automatic proximal femur bone contour extraction on clinical data

AAM/ASM by simultaneously optimizing both the global shape constraints and local image feature information. Experiments on clinical data sets verified the validity and performance of this approach. A further observation is that a traditional ICP/ASM alike method is used in our 2D/3D registration to fit the 3D model with the extracted contours. Due to the fact that a proper initial alignment of the 3D model is available, this works well in our approach. It can be expected that introducing the graphical model based Bayesian inference in the process to find the 2D/3D correspondence and the 3D model parameters may further improve its robustness and accuracy. R EFERENCES [1] Y. Chen, X. H. Ee, W. K. Leow, and T. S. Howe, “Automatic extraction of femur contours from hip x-ray images,”


in CVBIA’05, 2005, pp. 200–209. [2] R. de Luis-Garcia, M. Martin-Fernadez, J. I. Arribas, and C. Alberola-Lopez, “A fully automatic algorithm for contour detection of bones in hand radiographs using active contours,” in ICIP’03, Vol. 3, 2003, pp. 421–424. [3] M. G. Roberts, T. F. Cootes, and J. E. Adams, “Automatic segmentation of lumbar vertebrae on digitised radiographs using linked active apperance models,” in MIUA’06, 2006, pp. 120–124. [4] H. Gottschling, M. Roth, A. Schweikard, , and R. Burgkart, “Intraoperative, Fluoroscopy-based Planning for Complex Osteotomies of the Proximal Femur,” The International Journal of Medical Robotics & Computer Assisted Surgery, vol. 1, no. 3, pp. 33–38, 2005. [5] F. Bartolini, M. Carfagni, and L. Governi, “Model-based Extraction of Femoral Medulla Ducts from Radiographic Images,” Image and Vision Computing, vol. 22, pp. 173– 182, 2004. [6] T. P. Tian, Y. Chen, W. K. Leow, W. Hsu, T. S. Howe, and M. A. Png, “Computing neck-shaft angle of femur for x-ray fracture detection,” in CAIP’03, 2003, pp. 82–89. [7] S. Benameur, M. Mignotte, S. Parent, H. Labelle, W. Skalli, and J. A. de Guise, “3D/2D registration and segmentation of scoliotic vertebrae using statistical models,” Computerized Medical Imaging and Graphics, vol. 27(5), pp. 532–539, 2003. [8] H. Lamecker, T. H. Wenckebach, and H. C. Hege, “Atlasbased 3d-shape reconstruction from x-ray images,” in ICPR’06, Vol. I, 2006, pp. 371–374. [9] T. S. Tang and R. E. Ellis, “2d/3d deformable registration using a hybrid atlas,” in MICCAI’05, Vol. 3750, 2005, pp. 223–230. [10] G. Behiels, D. V. der Meulen, F. Maes, P. Suetens, and P. Dewaele, “Active shape model-based segmentation of digital x-ray images,” in MICCAI’99, 1999, pp. 128–137. [11] B. Howe, A. Gururajan, and L. R. Long, “Hierarchical segmentation of cervical and lumbar vertebrae using a customised generalized hough transform and extensions to active appearance models,” in Proc. IEEE 6th SSIAI, 2004, pp. 182–186. [12] G. Langs, P. Peloschek, and H. Bischof, “Determining position and fine shape detail in radiological anatomy,” in Pattern Recognition’03, LNCS 2781, 2003, pp. 532–539. [13] M. Seise, S. J. McKenna, I. W. Ricketts, and C. A. Wigderowitz, “Probabilistic segmentation of the knee joint from x-ray images,” in MIUA’06, 2006, pp. 110–114. [14] M. de Bruijne and M. Nielsen, “Image segmentation by shape particle filtering,” in ICPR’04, Vol. 3, 2003, pp. 722– 725. [15] B. Ma and R. E. Ellis, “Surface-based registration with a particle filter,” in MICCAI’04, 2004, pp. 566–573. [16] G. Zheng and L.-P. Nolte, “Surface reconstruction of bone from x-ray images and point distribution model incorporating a novel method for 2d-3d correspondence,” in CVPR’06, 2006, pp. 2237–2244. [17] C. Andrieu and A. Doucet, “Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC,” IEEE Trans. Signal Processing, vol. 47, no. 10, pp. 2667–2676, 1999. [18] M. Briers, A. Doucet, and S. Maskell, “Smoothing algorithms for state-space models,” Technical Report CUED/FINFENG/TR.498. Cambridge University Engineering Department, 2004. [19] J. Yedidia, W. Freeman, and Y. Weiss, “Understanding belief propagation and its generalizations, ISBN: 1-55860811-7,” Exploring artificial intelligence in the new millennium, pp. 239–269, 2003. [20] A. P. Young, “Spin glasses: A computational challenge for the 21st century.” [Online]. Available: citeseer.ist.psu.edu/young01spin.html


53

[21] G. Parisi, “A backtracking survey propagation algorithm for K-satisfiability,” arXiv:cond-mat/0308510v1, 2003. [22] T. S. Caetano, T. Caeli, and D. A. C. Barone, “An optimal probabilistic graphical model for point set matching,” Technical Report TR 04-03, University of Alberta, Edmonton, Alberta Canada, 2004, unpublished. [23] J. Larocque, J. P. Reilly, and W. Ng, “PArticle filters for tracking an unknown number of sources,” IEEE Trans. Signal Processing, vol. 50, no. 12, pp. 2926–2937, 2002. [24] R. J. McEliece, D. J. C. MacKay, and J. F. Cheng, “Turbo decoding as an instance of Pearl’s beliefpropagation algorithm,” IEEE Journal on Selected Areas in Communications, vol. 16, no. 2, pp. 140–152, 1998. [25] T. J. Richardson and R. L. Urbanke, “The capacity of lowdensity parity-check codes under message-passing decoding,” IEEE Trans. on Information Theory, vol. 47, no. 2, pp. 599–618, 2001. [26] T. Minka, “A family of algorithms for approximate Bayesian inference,” PhD Thesis, 2001. [27] J. M. Pena, J. A. Lozano, and P. Larranaga, “Unsupervised learning of Bayesian networks via estimation of distribution algorithms: an application to gene expression data clustering,” International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems, vol. 12, pp. 63–82, 2003. [28] M. Isard, “A smoothing filter for Condensation,” in ECCV’98, 1998, pp. 767–781. [29] L. Sigal, S. Bhatia, S. Roth, M. J. Black, and M. Isard, “Tracking loose-limbed people,” in CVPR’04, 2004, pp. 421–428. [30] Y. Wu, G. Hua, and T. Yu, “Tracking articulated body by dynamic markov network,” in ICCV’03, 2003, pp. 1094– 1101. [31] E. B. Sudderth, M. I. Mandel, W. T. Freeman, and A. S. Willsky, “Visual hand tracing using nonparametric belief propagation,” in CVPRW’04, Vol. 12, 2004, pp. 189–197. [32] X. Dong and G. Zheng, “Fully automatic determination of morphological parameters of proximal femur from calibrated fluoroscopic images thourgh particle filtering,” in ICIAR’06, 2006, pp. 182–186. [33] G. Hua and Y. Wu, “Variational maximum a posteriori by annealed mean field analysis,” IEEE Trans. on Pattern Analysis and Machine Intelligence, vol. 27, no. 11, pp. 1747–1761, 11 2005. [34] M. Isard, “Pampas: Real-valued graphical models for computer vision,” in CVPR’03, 2003, pp. 613–620. [35] M. W. Lee and I. Cohen, “Human upper body pose estimation in static images,” in ECCV’04, 2004, pp. 126– 138. [36] J. Coughlan and S. Ferreira, “Finding deformable shapes using loopy belief propagation,” in ECCV’02, 2002, pp. 453–468. [37] A. Rangarajan, J. Coughlan, and A. L. Yuille, “A bayesian network framework for relational shape matching,” in ICCV’03, 2003, pp. 671–678. [38] V. C. Raykar and R. Duraiswami, “Fast optimal bandwidth selection for kernel density estimation,” in Proceedings of the 6th SIAM International Conference on Data Mining, 2006, pp. 524–528. [39] J. MacCormick and M. Isard, “Partitioned Sampling, Articulated Objects, and Interface-Quality Hand Tracking,” in ECCV’00, 2000, pp. 3–19.

54

Xiao Dong is a PhD candidate in the Institute for Surgical Technology and Biomechanics, MEM Research Center, University of Bern, Switzerland. He received his Masters in 1998 from Beijing University of Posts and Telecommunication (BUPT), China. From 2000, he worked in the Electrical Engineering Department, Eindhoven University of Technology as a research assistant and received his PDEng certificate in 2002. He is currently working on the statistical model based automatic identification, segmentation, and reconstruction of anatomical structures from calibrated fluoroscopic images.

Miguel A. Gonzalez Ballester obtained a PhD from the University of Oxford (Robotics Research Group) in 1999, in the area of medical image analysis for the study of brain morphology and its link to the aetiology of schizophrenia and multiple sclerosis. He then worked for two years in Thshiba Medical Systems in Otawara, Japan, doing research for their MRI deivision. In Dec. 2001 he obtained a faculty position at the Epidaure (currently renamed Asclepios) Research Group in INRIA at Sophia Antipolis, France, where he worked on image registration and statistical methods for medical image analysis, until he moved to Switzerland in mid 2004. Dr. González is currently leading the Surgical Technology Division of the Institute for Surgical Technology and Biomechanics, University of Bern, as well as heading the Medical Image Analysis Group.

Guoyan Zheng obtained his Ph.D. in biomedical engineering from University of Bern in 2002. He is currently a Research Scientist and the head of the Surgical Navigation Group of the same university at the Institute for Surgical Technology and Biomechanics, MEM Research Center. His research interests include fluoroscopy image calibration and processing, statistical model based segmentation and reconstruction, 2D/3D intensityand feature-based registration, endoscopy-based tracking, and computer-assisted interventions. His works are partially supported by the Swiss National Center for Competence in Research “Computer Aided and Image Guided Medical Interventions” (http://co-me/) located at the ETHZ and the Swiss Innovation Promotion Agency CTI (http://www.kti-cti.ch). He is a member of the IEEE, the CAOS-International, and the MICCAI society.

