In this paper an algorithm to carry out the automatic segmentation of bone structures in 3D CT images
has been implemented. Automatic segmentation of bone structures is of special interest for radiologists and
surgeons to analyze bone diseases or to plan some surgical interventions. This task is very complicated as bones
usually present intensities overlapping with those of surrounding tissues. This overlapping is mainly due to the
composition of bones and to the presence of some diseases such as Osteoarthritis, Osteoporosis, etc. Moreover,
segmentation of bone structures is a very time-consuming task due to the 3D essence of the bones. Usually,
this segmentation is implemented manually or with algorithms using simple techniques such as thresholding and
thus providing bad results. In this paper gray information and 3D statistical information have been combined
to be used as input to a continuous max-
ow algorithm. Twenty CT images have been tested and di erent
coe cients have been computed to assess the performance of our implementation. Dice and Sensitivity values
above 0.91 and 0.97 respectively were obtained. A comparison with Level Sets and thresholding techniques has
been carried out and our results outperformed them in terms of accuracy.Ministerio de ciencia e innovación TEC2010-21619-C04-02Junta de Andalucía P11-TIC-772

Acha Piñero, Begoña

Pérez Carrasco, José Antonio

Serrano Gotarredona, María del Carmen

English

idUS. Depósito de Investigación Universidad de Sevilla

Segmentation of Bone Structures in 3D CT Images Based onContinuous Max-flow OptimizationJ. A. Pe´rez-Carrascoa, B. Achaa and C. SerranoaaUniversity of Seville Signal and Communications Dpt. Camino de los Descubrimientos, s/n.41092,Seville,SpainABSTRACTIn this paper an algorithm to carry out the automatic segmentation of bone structures in 3D CT imageshas been implemented. Automatic segmentation of bone structures is of special interest for radiologists andsurgeons to analyze bone diseases or to plan some surgical interventions. This task is very complicated as bonesusually present intensities overlapping with those of surrounding tissues. This overlapping is mainly due to thecomposition of bones and to the presence of some diseases such as Osteoarthritis, Osteoporosis, etc. Moreover,segmentation of bone structures is a very time-consuming task due to the 3D essence of the bones. Usually,this segmentation is implemented manually or with algorithms using simple techniques such as thresholding andthus providing bad results. In this paper gray information and 3D statistical information have been combinedto be used as input to a continuous max-flow algorithm. Twenty CT images have been tested and differentcoefficients have been computed to assess the performance of our implementation. Dice and Sensitivity valuesabove 0.91 and 0.97 respectively were obtained. A comparison with Level Sets and thresholding techniques hasbeen carried out and our results outperformed them in terms of accuracy.Keywords: Segmentation, CT, Max-flow, Optimization1. INTRODUCTIONAutomatic and accurate segmentation of bones is of special interest for radiologists and surgeons in order toanalyze and locate many kinds of fractures, to diagnose some bone diseases (such as Osteoarthritis, RheumatoidArthritis, Osteoporosis) and in plastic surgical planning. The segmentation of such structures is difficult becausethe presence of some diseases usually provokes a reduction in bone density in different areas. Besides, the differentparts of the bones (mainly periosteum, compact (hard) bone, cancellous (spongy) bone and bone marrow) presentwide differences in densities and thus, different Hounsfield values.1 This implies an overlapping of these valueswith other tissue types such as muscle, fat or some organs.In Fig. 1 several bone structures can be seen. It can be seen in the figure that the presence of a widerange of different Hounsfield values in these structures make the segmentation problem more coplicated (notethe different Hounsfiel values in cortical and cancellous bone).A second difficulty in the segmentation of bone structures is their 3D essence. Thus, fast and automaticalgorithms are required in order to obtain the segmentation with reduced computational times.Many algorithms have been published regarding to bone segmentation and, according to Wang et al.,2they can be classified into four categories: intensity-based,3 edge-based,4 region-based,5 and model-based.4,6In some works the segmentation technique is selected according to the number of bones or type of bone tobe segmented.7–9 This lack of generalization is a problem, as radiologists do not want to consider differentsegmentation methods according to the regions to be analyzed. Thus, a more general method more suitablefor a wide number of bone structures would be desirable. When considering more general methods, level sets,active contours and graph cuts implementations are preferred.10–12 However, level set methods can drive toFurther author information: (Send correspondence to J. A. Pe´rez-Carrasco)J. A. Pe´rez-Carrasco: E-mail: jperez2@us.es, Telephone: +34 954 48 73 33B. Acha: E-mail: bacha@us.es, Telephone: +34 954 48 73 33C. Serrano: E-mail: cserrano@us.es, Telephone: +34 954 48 73 33Medical Imaging 2015: Image Processing, edited by Sébastien Ourselin, Martin A. Styner, Proc. of SPIE Vol. 9413, 94133Y · © 2015 SPIE · CCC code: 1605-7422/15/$18 · doi: 10.1117/12.2082139Proc. of SPIE Vol. 9413  94133Y-1Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 1/18/2019Terms of Use: https://www.spiedigitallibrary.org/terms-of-useCompact BoneCancellous BoneFigure 1. Zoomed abdominal section of a CT volume. Note differences in compact bone with cancellous bone in thedifferent bone structures.local optima of the minimization energy function and suffer from high sensitivity to initialization. Graph cuttechniques have the advantage that they guarantee the global optima (if bipartition is implemented and a Pottsmodel is adopted) in nearly real time.In the algorithm presented in this paper, a methodology based on continuous convex relaxation tech-niques13,14 is followed. We combine gray level information, gradient values and statistical information computedusing histogram distances and provide this combined information as input to a convex relaxation algorithm inorder to get accurate and fast segmentation of bones in 3D volumes. Continuous convex relaxation techniquesshare the advantages of both active curves and graph-cuts14 and they outperform graph cuts in regard to speedand accuracy.142. METHODThe methodology implemented in the work presented in this paper consists in three steps. The first step isa bone enhancement operation. Hounsfield values in CT volumes are usually between -2000 and 3000. Thebones usually have values from 700 (cancellous bone) to 3000 (dense bone) approximately. In order to get ahigher contrast, a thresholding operation followed by a normalization operation has been applied to the differentimages used throughout our experiments. The resulting normalized image is denoted as Inorm. The threshold,maximum and minimum values used during the two operations are common in all the volumes that have beenused in order to not lose generalization.The second stage of our algorithm is the computation of Histogram Distance Images (HDI ). 10% of thevolumes used throughout our experiments have been used to obtain a model histogram of bone structures. Therest 90% of the volumes were used to test our algorithm and for comparison purposes. In this dataset, for eachvoxel (i, j, k), a local histogram is computed. Different 3D sizes of this local neighborhood were tested and alocal neighborhood of 7x7x3 provided the best results. Subsequently, for each local histogram, a distance17 tothe histogram model is calculated and this distance is assigned to voxel (i, j, k) in a new Histogram DistanceImage (HDI ). HDI will have values close to 0 in regions where bone tissue appears and close to 1 when nobones are encountered. Finally, in order to combine this 3D statistic information with the gray level informationInorm, we create a cost image term to be minimized (CIT ) as:CIT = ((1− Inorm) +HDI)/2 (1)Eq. 1 is the cost to be minimized using the continuous max-flow algorithm described in the next subsection.The first term of the right expression in Eq. 1 uses intensity information. Pixels corresponding to bone (mainlycortical bone) will have values close to ‘1’. Thus, we use 1− Inorm instead of Inorm directly because we wantlow values in bone structures. Something similar occurs with the second term in the right expression in Eq. 1.The distance to the model histogram in pixels corresponding to bones should be close to 0 whereas it shouldbe close to 1 in pixels belonging to other kind of tissues. In Fig. 2 these computed images are shown for a CTProc. of SPIE Vol. 9413  94133Y-2Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 1/18/2019Terms of Use: https://www.spiedigitallibrary.org/terms-of-usea)e)b)(1)Figure 2. Computed images used by the algorithm. The first image (a) is the initial image without thresholding. Thesecond image (b) corresponds to the thresholded and normalized image. The third image (c) corresponds to the distanceimage HDI. Finally, the last image (d) shows the new computed image CIT.slice. The first image (a) is the initial image without thresholding. The second image (b) corresponds to thethresholded and normalized image. The third image (c) corresponds to the distance image HDI. Finally, thelast image (d) corresponds to CIT image which will be the input to the continuous max-flow algorithm. It canbe seen that bone structures are darker in CIT image, which means that gray level information and distancebetween histograms information have been properly combined.Finally, in the last and third stage, we use the CIT image above computed, as input to the continuousmax-flow optimization algorithm proposed by J. Yuan et al.13 In the work by Yuan et al., citeYuan20102217,Yuan et al. proved that the max-flow problem, which is formulated by maximizing the total flow from thesource, is equivalent to a continuous s-t min-cut problem as follows:minu(x)[0,1]∫Ω(u)Csdx+∫Ω(1− u)Ctdx+∫ΩC(x)|∇u|dx (2)u(x) is the labeling function and indicates if point x belongs to the region to be segmented. If the minimizationproblem is well defined, the cost function Cs should take low values inside the bones and high values outsidethem. Similarly, Ct should take low values outside the bones and high values inside them. The most right termof Eq. (2) is the regularization term and C(x) is a penalty function. |∇u| is the absolute gradient of the labelingfunction u(x), thus indicating the boundary of the segmented region. Cs and Ct are called regional terms. Inthe algorithm proposed here, the regional term Cs, provided as input to the continuous max-flow algorithm, isProc. of SPIE Vol. 9413  94133Y-3Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 1/18/2019Terms of Use: https://www.spiedigitallibrary.org/terms-of-usethe CIT image computed through Eq. (1). Thus, in the algorithm proposed, the terms Cs and Ct (Eq. 2) arecomputed as follows:Cs = CIT (3)Ct = 1− CIT (4)Note that with Eq. (4) we force Ct to be low outside the bones and high inside them, as required by theminimization in (2). In addition, the algorithm penalizes an excessive area of the surface limiting the segmentedvolume. This penalization term is the most right term in Eq. (2). This penalization is modeled by a functionC(x) depending on the gradient of the CIT image along the surface, so that if the gradient is high along thesurface, a large area of the surface is not penalized. The penalization function is computed as follows:C(x) =b1 + a|∇CIT (x)| (5)where parameters a and b control the importance of the gradient in the penalization function. In our imple-mentation a and b values have been obtained empirically and their values are 10 and 0.2, respectively. Notethat at the surface, |∇u| will have high values whereas C(x) will have low values if the gradient of CIT (x) ishigh. Thus, the most right term in Eq. 2 will not penalize a large area if a high border is encountered.Finally, in order to smooth out the resulting segmented image, and to get a more accurate segmentation, amorphological operation of opening is performed.It is important to note that the proposed algorithm uses the gradient information that is very high inthe boundaries of cortical bone structures and histogram distance information, which is very useful mainly tosegment cancellous bone structures, which present lower gray values and overlap with other tissue types suchas muscle, fat or some organs.3. EXPERIMENTAL RESULTSTo assess the algorithm proposed, twenty images corresponding to ten different patients have been used. Manualsegmentations of all the slices were provided by an expert. Three different methods have been compared usingthe same set of images. Table 1 shows the performance parameters provided by each of the techniques. Theseparameters are the Jaccard, Dice, Sensitivity (S ) and Positive Predictive Value (PPV ) parameters. Both Jaccardand DICE coefficients measure the set agreement in terms of false positive, false negative, true negative andtrue positive counts.Table 1. Comparison between the three different algorithms.METHOD Computational Time PPV Sensitivity Specificity Dice JaccardTHRESHOLDING 0,017 s 0,8 0,802 0,993 0,7857 0,656DRLSE 3402 s 0,54 0,938 0,889 0,6031 0,511thresholded HDI 255 s 0.872 0.9570 0.9946 0.9106 0.8397The first method was a technique based on thresholding. A common threshold was computed for all theimages and the threshold providing the best results was 1200. As it has been stated before, the low resultsprovided by this method are because bone structures do not have a pixel-wise constant object distribution,which means that they are composed by different structures with different Hounsfield values.The second method employed has been the Distance Regularized Level Set Evolution (DRLSE ) published inRef. 18. In the experiments shown in Ref. 18, low computational times were obtained when small images witha reduced number of regions were employed. However, when bigger images are considered and many differentbone structures are considered, performance values were considerably low and excessively high computationaltimes were required. In all these cases the algorithm used to have convergence problems. Table 1 shows theresults when using Inorm image as input and using parameters λ = 20 , timestep ∆t = 5,  = 1.5 and α = 1 .Proc. of SPIE Vol. 9413  94133Y-4Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 1/18/2019Terms of Use: https://www.spiedigitallibrary.org/terms-of-use(tr-0Ar"YFigure 3. Results obtained for four different images using our algorithm. Upper row shows the original images. Lower rowshows the corresponding groundtruth images. Results obtained with the DRLSE algorithm are shown in red. Resultsprovided by our algorithm are shown in blue.The third method corresponds to our algorithm and it has been labeled as thresholded HDI. Only pixelsabove a threshold value were considered to compute the HDI image. Using this thresholding operation areduction of 80% in the computational time was obtained when compared with the same technique withoutthresholding.Fig. 3 shows the results obtained for four different images using our algorithm and DRLSE. Note how theDRLSE algorithm has problems when some bone structures have inner regions that do not correspond to bones.The segmentation step corresponding to the continuous max-flow algorithm offered an average computationalcost of 0.5s per 512x512 slice, which is much smaller than values provided by other minimization strategies.The creation of the HDI image was the most time consuming task (see Table 1). Matlab 7.12 was used in allthe experiments.4. CONCLUSIONSThe segmentation of bone structures is a complicated task because they present intensities overlapping withthose of surrounding tissues. In this paper an automatic algorithm, which combines gray level information andstatistical information extracted from local histograms in a 3D neighborhood of the voxels under consideration,is computed for segmentation of bone structures. The gray information and the statistical information arecombined and used as input to a continuous max-flow optimization algorithm. Several coefficients have beencomputed to measure the algorithm performance. Results obtained are over 0.9, which indicates neither theFP rate nor the FN rate are high. High sensitivity values above 0.97 were obtained. A further validation,with a higher number cases would be desirable. Note that the segmentation stage using a continuous max-flowimplementation is faster than the majority of algorithms in the literature. Future implementations will speedup the creation of the input image CIT, which is the main bottleneck in our implementation.The authors declare that the work is not being, or not has been, submitted for publication or presentationelsewhere.ACKNOWLEDGMENTSThis research has been supported by TEC2010-21619-C04-02 and P11-TIC-7727.Proc. of SPIE Vol. 9413  94133Y-5Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 1/18/2019Terms of Use: https://www.spiedigitallibrary.org/terms-of-useREFERENCES[1] Feeman, T. G., [The mathematics of medical imaging: A beginner’s guide], Springer undergraduate texts inmathematics and technology. Springer (2010).[2] Wang, L. G., et al., “Validation of bone segmentation and improved 3-D registration using contour coherencyin CT data,” IEEE Transactions on Medical Imaging, 25(3), 324–334, (2006).[3] Kang, Y., et al., “A new accurate and precise 3-D segmentation method for skeletal structures in volumetricCT data,” IEEE Transactions on Medical Imaging, 22(5), 586–598, (2003).[4] Pardo, X. M., et al., “A snake for CT image segmentation integrating region and edge information,” Imageand Vision Computing, 19(7), 461–475, (2001).[5] Beveridge, J. R., et al., “Segmenting images using localized histograms and region merging,”InternationalJournal of Computer Vision, 2(3), 311–347, (1989).[6] Sebastian, T. B., et al., “Segmentation of carpal bones from CT images using skeletally coupled deformablemodels,” Medical Image Analysis, 7(1), 2–45, (2003).[7] Cheng, Y., et al., “Automatic segmentation technique for acetabulum and femoral head in CT im-ages,”Pattern Recognition, 46(11), 2969–2984, (2013).[8] Cervinka, T., Hyttinen, J., Sievnen, H., “Accurate cortical bone detection in peripheral quantitative com-puted tomography images,” in IFMBE Proceedings, 289–292, (2014).[9] Aslan, M. S., et al., “A novel 3D segmentation of vertebral bones from volumetric CT images using graphcuts,” Lecture Notes in Computer Science, 5876 LNCS, 519–528, (2009).[10] Shadid, W., Willis, A., “Bone fragment segmentation from 3D CT imagery using the Probabilistic Water-shed Transform,” in Conference Proceedings - IEEE SOUTHEASTCON, (2013).[11] Calder, J., et al., “A variational approach to bone segmentation in CT images,” in Progress in BiomedicalOptics and Imaging, SPIE (2011).[12] Kratky, J., “Three-dimensional segmentation of bones from CT and MRI using fast level sets,” in Progressin Biomedical Optics and Imaging, SPIE, (2008).[13] Yuan, J., Bae, E., Tai, X. C., “A study on continuous max-flow and min-cut approaches,” in Proceedingsof CVPR, 2217–2224, (2010).[14] Punithakumar, K., Yuan, J., “A Convex Max-Flow Approach to Distribution-Based Figure-Ground Sepa-ration,” SIAM Journal of Imaging Sciences, 5(4),1333–1354, (2012).[15] Chan, T. F., et al., “Algorithms for finding global minimizers of image segmentation and denoising models,”SIAM Journal on Applied Mathematics, 66(5), 1632–648, (2006).[16] Bresson, X., et al., “Fast global minimization of the active contour/snake model,” Journal of MathematicalImaging and Vision, 28(2), 151–167, (2007).[17] Rubner, Y., et al., “The Earth Mover’s Distance as a Metric for Image Retrieval,” International Journalof Computer Vision, 40(2), 99–121, (2000).[18] Li, C., Xu, C., Gui, C., Fox, M. D., “Distance Regularized Level Set Evolution and its Application toImage Segmentation,” IEEE Trans. Image Processing, 19(12), 3243–3254, (2010).Proc. of SPIE Vol. 9413  94133Y-6Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 1/18/2019Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Segmentation of bone structures in 3D CT images based on continuous max- ow optimization

https://idus.us.es/bitstream/handle/11441/81749/94133Y.pdf?sequence=1&isAllowed=y

Segmentation of bone structures in 3D CT images based on continuous max- ow optimization

Abstract

Similar works

Full text

Available Versions

idUS. Depósito de Investigación Universidad de Sevilla