SEGMENTATION OF TRABECULAR BONES ... - Semantic Scholar

Report 3 Downloads 96 Views
SEGMENTATION OF TRABECULAR BONES FROM VERTEBRAL BODIES IN VOLUMETRIC CT SPINE IMAGES Melih S. Aslan, Asem Ali, Ben Arnold∗ , Rachid Fahmi, Aly A. Farag, and Ping Xiang ∗ University of Louisville Computer Vision and Image Processing Laboratory Louisville, KY, 40299, USA



Image Analysis, Inc. 1380 Burkesville St. Columbia, KY, 42728, USA

ABSTRACT We present a 3D segmentation technique of trabecular (cancellous) bones in CT images of Vertebral bodies (VBs). In order to be used for Bone Mineral Density (BMD) measurements, the cortical and trabecular bones are subsequently segmented using graph cuts method and local volume growing methods separately. In the final step, we measure our segmentation accuracy for each method. Validity was analyzed using ground truths of data sets and the European Spine Phantom (ESP). Preliminary results are very encouraging and a reproducibility of the results was achieved for 16 data sets. The average segmentation error is below 2.0% for both methods. Index Terms— Spine column, trabecular bone, 3D Segmentation, Bone Mineral Density, segmentation accuracy. 1. INTRODUCTION In this paper, we are primarily interested in volumetric computed tomography (CT) images of the vertebral bones of spine column with a particular focus on the lumbar spine. The primary goal of the proposed work is in the field of spine densitometry where bone mineral density (BMD) measurements are restricted to the vertebral bodies (see Fig. 1 for regions of spine bone). Various approaches have been introduced to tackle the 3D segmentation of skeletal structures in general and of vertebral bodies in particular for the anatomical definition of a VB). For instance, Kang et al. [1] proposed a 3D segmentation method for skeletal structures from CT data. Their method is a multi-step method that starts with a three dimensional region growing step using local adaptive thresholds followed by a closing of boundary discontinuities and then an anatomically-oriented boundary adjustment. Applications of this method to various anatomical bony structures are presented and the segmentation accuracy was determined using the European Spine Phantom (ESP) [2]. Later, Mastmeyer et al. [3] presented a hierarchical segmentation approach for the lumbar spine in order to measure bone mineral density. This approach starts by separating the vertebrae from each other followed by a two step segmentation using a deformable

Fig. 1. Anatomy of a human vertebra (The image is adopted from [4]).

mesh followed by adaptive volume growing operations. The authors conducted a performance analysis using two phantoms: a digital phantom based on an expert manual segmentation and the ESP. They also reported that their algorithm can be used to analyze three vertebrae in less than 10min. This timing is far from the real time required for clinical applications but it is a huge improvement compared to the timing of 1 − 2h reported in [5]. Recently, in the context of evaluating the Ankylosing Spondylitis, Tan et al. [6, 7] presented a technique to segment whole vertebrae with their syndesmophytes using a 3D multi-scale cascade of successive level sets. The seed placement was done manually and results were validated using synthetic and real data. The implementation of this technique uses the open library ITK functions [8]. Other techniques have been developed to segment skeletal structures and can be found for instance in [9, 10] and the references therein. The VB consists of trabecular and cortical bones. The main objective of our algorithm is the segmentation of the trabecular region of the vertebral body using a graph cuts seg-

mentation method proposed by Asem and Farag [11] and local volume growing algorithm proposed by Kang et al. [1]. The segmentation results of trabecular bones may be used to determine various VOIs for bone mineral density measurements (see for instance [3, 12, 13]). Section 2.1 discusses the background of Graph Cuts method [11]. Section 2.2 discusses the local adaptive volume growing algorithm [1]. 2. MATERIALS AND METHODS In this section, we briefly review the mathematical background of the techniques used in this paper. We first describe the basic idea of graph cuts. In the second part, we describe briefly local adaptive volume growing method. 2.1. Graph Cuts Segmentation Framework To segment a trabecular bone, we initially labeled the volume based on its gray level probabilistic model. Then we create a weighted undirected graph with vertices corresponding to the set of volume voxels P, and a set of edges connecting these vertices. Each edge is assigned a nonnegative weight. The graph also contains two special terminal vertices s (source) “trabecular bone”, and t (sink) “cortex”. Consider a neighborhood system in P, which is represented by a set N of all unordered pairs {p, q} of neighboring voxels in P. Let L the set of labels {“0”, “1”}, correspond to trabecular and cortex regions respectively. Labeling is a mapping from P to L, and we denote the set of labeling by f = {f1 , . . . , fp , . . . , f|P| }. In other words, the label fp , which is assigned to the voxel p ∈ P, segments it to trabecular or cortex region. Now our goal is to find the optimal segmentation, best labeling f , by minimizing the following energy function: X X E(f ) = Dp (fp ) + V (fp , fq ), (1) p∈P

{p,q}∈N

C+

P (Ip |fp ) =

r=1

µ γ∗ =

K−

¶ K2 fneq (f ) . K −1

(3)

where K = 2 is the number of classes in the volume and fneq (f ) denotes the relative frequency of the not equal labels in the voxel pairs. To segment a volume, instead of independently segmenting each 2D slice of the volume, we segment the 3D trabecular bone using a 3D graph where each vertex in this graph represents a voxel in the trabecular volume. Then we define the weight of each edge as shown in table 2.1. After that, we get the optimal segmentation surface between the trabecular and its background by finding the minimum cost cut on this graph. The minimum cost cut is computed exactly in polynomial time for two terminal graph cuts with positive edges weights via s/t Min-Cut/Max-Flow algorithm [15].

Edge {p, q} {s, p} {p, t}

Weight γ 0 −ln[P (Ip | “1”)] −ln[P (Ip | “0”)]

for fp 6= fq fp = fq p∈P p∈P

2.2. A Locally Adaptive Volume Growing Framework

where Dp (fp ), measures how much assigning a label fp to voxel p disagrees with the voxel intensity, Ip . Dp (fp ) = −ln P (Ip | fp ) is formulated to represent the regional properties of segments. The second term is the pairwise interaction model which represents the penalty for the discontinuity between voxels p and q. To initially label the trabecular volume and to compute the data penalty term Dp (fp ), we use the modified EM [14] to approximate the gray level marginal density of each class fp , trabecular and cortex region, using a LCG with Cf+p positive and Cf−p negative components as follows: fp X

positive weight in class fp and wf−p ,l means the lth negative weight in class fp . This weights have a restriction + − PCfp PCfp + − r=1 wfp ,r − l=1 wfp ,l = 1. The simplest model of spatial interaction is the Markov Gibbs random field (MGRF) with the nearest 6-neighborhood. Therefore, for this specific model the Gibbs potential can be obtained analytically using our maximum likelihood estimator (MLE) for a generic MGRF in [11]. So, the resulting approximate MLE of γ is:

Three dimensional region growing or simply volume growing algorithms use global criteria to segment components of an image [16]. However, such algorithms do not yield satisfactory results for the spine CT images due to spatial resolution, noise, as well as to bone deformations and diseases. In this work, we use a locally adaptive criterion [1]. This criterion is based on the mean intensity value and its standard deviation in the 26-neighborhood of the considered voxel as follows: ( if I(v) ≥ µ − ασ, label v as cortical if I(v) < µ − ασ, label v as trabecular

(4)

C−

wf+p ,r ϕ(Ip |θf+p ,r ) −

fp X

wf−p ,l ϕ(Ip |θf−p ,l ),

l=1

(2) where ϕ(.|θ) is a Gaussian density with parameter θ ≡ (µ, σ 2 ) with mean µ and variance σ 2 . wf+p ,r means the rth

The procedure of this algorithm is to check every voxel on the outer surfaces and if its gray value is greater than a local threshold then use it to start a local volume growing. The separated regions define the VOIs that will be used for integral and trabecular density measurements respectively.

Table 1. Accuracy of segmentation of trabecular regions on 16 data sets. Average volume 512x512x15.

(a)

(b)

(c)

Min. error, % Max. error, % Mean error, % Stand. dev.,% Average time, sec

Algorithm Graph Cuts [11] Volume Growing [1] 0.05 0.03 5.01 2.56 1.07 0.67 1.21 0.69 17.6 19.2

Fig. 2. (a) Integral segmented VB including the trabecular and cortical bones. (b) Outline of trabecular bone separated from cortical shell. (c) A sample 3D representation of trabecular (inner region) and cortex (outer region) bone. 3. EXPERIMENTS AND DISCUSSION The main goal of our experiment is the segmentation trabecular region of vertebral bone as shown in Fig 2. We use two different methods separately. The first method is Graph Cuts based segmentation framework proposed by Asem and Farag [11] and the second method is volume growing model proposed by Kang et al. [1]. At the end of our experiment, we get two segmented trabecular bone images of each VB to be compared with its ground truth. To evaluate the results we calculate the percentage segmentation error as follows:

error% =

100 ∗ N umber of missclassif ied voxels (5) N umber of trabecular bone voxels

Preliminary results are very encouraging and a reproducibility of the results was achieved for 16 data sets, which we have their ground truths (expert segmentation). The statistical analysis of data sets are shown in Table 1. In this table the results of both graph cut segmentation [11] and local adaptive volume growing method [1] are shown.The motivation behind our segmentation of trabecular bone is to exclude errors as far as possible. All algorithms are un on a PC 3Ghz AMD Athlon 64 X2 Dual, and 3GB RAM. All implementations are in C++. 4. VALIDATION To assess the robustness of the proposed approach, in addition to ground truths of data sets, the European Spine Phantom (ESP), which is an accepted standard for quality control [2] in bone densitometry, was used to validate the segmentation algorithms. Using Siemens Sensation 64 generation, ESP was scanned at 120kV and 1mm slice thickness. Fig. 3 shows some slices of the vertebral volume of ESP. The error of the trabecular bone segmentation is 1.05% and 1.86% for graph

Fig. 3. Some CT images from the ESP data set. cuts method and adaptive volume growing methods, respectively. Fig. 4 shows 3D segmentation results. 5. CONCLUSION We have presented a 3D segmentation framework of trabecular bones in CT images. In order to be used for BMD measurements, trabecular bone is separated using graph cuts and volume growing approaches, respectively. The main objective of this work is its application for an ongoing BMD measurement project in both the thoracic and lumbar spines. Preliminary results are very encouraging and a reproducibility of the results was achieved for 16 data sets. The segmentation error is below 2% for both methods. 6. ACKNOWLEDGEMENTS This research has been supported by Image Analysis, Inc in Columbia, KY. Special thanks go to Dr. Dongqing Chen who spent his valuable time for technical discussion in this experiment. 7. REFERENCES [1] Y. Kang and K. Engelke and W. A. Kalender, New accurate and precise 3D segmentation method for skeletal structures in volumetric CT data IEEE Transaction on Medical Imaging (TMI). vol. 22, no. 5, pp. 586-598, 2003.

bones from CT images using skeletally coupled deformable models Medical Image Analysis vol.7, no.1, pp. 21-45, 2003. [10] L. C. Perelman and J. Paradis and E. Barrett, R. A. Zoroofi and Y. Sato and T. Sasama and T. Nishii and N. Sugano and K. Yonenobu and H. Yoshikawa and T. Ochi and S. Tamura Automated segmentation of acetabulum and femoral head from 3-D CT images. IEEE Trans Inf Technol Biomed. vol. 7, no. 4, pp. 329-343, 2003.

(a)

(b)

Fig. 4. European Spine Phantom (ESP) segmentation results (dark area shows misclassified voxels) (a) The graph cut segmentation error is 1.05%, and (b) The volume growing segmentation error is 1.86 [2] W. A. Kalender and D. Felsenberg and H. Genant and M. Fischer and J. Dequeker and J. Reeve, The European Spine Phantom - a tool for standardization and quality control in spinal bone measurements by DXA and QCT European J. Radiology vol. 20, pp. 83-92, 1995. [3] A. Mastmeyer and K. Engelke and C. Fuchs and W. A. Kalender, A hierarchical 3D segmentation method and the definition of vertebral body coordinate systems for QCT of the lumbar spine Medical Image Analysis. vol. 10, no. 4, pp. 560-577, 2006. [4] www.back.com. [5] J. Kaminsky and P. Klinge and M. Bokemeyer and W. Luedemann and M Samii, Specially adapted interactive tools for an improved 3D-segmentation of the spine Computerized Medical Imaging and Graphics vol. 28, no. 3, pp. 119-127, 2004. [6] S. Tan and J. Yao and M. M. Ward and L. Yao and R. M. Summers, Computer aided evaluation of Ankylosing Spondylitis IEEE International Symposium on Biomedical Imaging (ISBI’06) pp. 339-342, 2006. [7] S. Tan and J. Yao and M. M. Ward and L. Yao and R. M. Summers, Computer Aided Evaluation of Ankylosing Spondylitis Using High-Resolution CT IEEE Transaction on Medical Imaging (TMI). vol. 27, no. 9, pp. 1252-1267, 2008. [8] L. Ibanez and W. Schroeder and J. Cates, The ITK software guide Kitware Inc., 2003. [9] Thomas B. Sebastiana and Hseyin Teka and Joseph J. Criscob and Benjamin B. Kimia, Segmentation of carpal

[11] A. M. Ali and A. A. Farag, Automatic Lung Segmentation of Volumetric Low-Dose CT Scans Using Graph Cuts Proc. of International Symposium on Visual Computing (ISVC’08) pp. 258-267, 2006. [12] W. A. Kalender and E. Klotz and C. Suess, Vertebral bone mineral analysis: an integrated approach with CT Radiology. vol. 164, no. 2, pp. 419-423, 1987. [13] P. Steiger and J. E. Block and S. Steiger and A. F. Heuck and A. Friedlander and B. Ettinger and S. T. Harris and C. C. Gluer and H. K. Genant, Spinal bone mineral density measured with quantitative CT: effect of region of interest, vertebral level, and technique Radiology vol. 175, no. 2, pp. 537-543, 1990. [14] A.A. Farag, A. El-Baz, G. Gimelfarb, Density estimation using modified expectation maximization for a linear combination of gaussians. Proc. of International Conference on Image Processing (ICIP’04) vol. 3, pp. 18711874, 2004. [15] Y. Boykol, V. Kolmogorov, An experiment comparison of min-cut/max-flowalgorithms for energy minimization in vision IEEE Transaction on Pattern Analysis and Machine Intelligence vol. 26, pp. 1124-1137, 2004. [16] R. Adams and L. Bischof, Seeded Region Growing IEEE Trans. Pattern Anal. Mach. Intell. vol. 16, no. 6, pp. 641-647, 1994.