ORIGINAL ARTICLES AND REVIEWS
Technique for bone volume measurement from human femur head samples by classification of micro-CT image histograms
Franco MarinozziI; Fabiano BiniI; Andrea MarinozziII; Francesca ZuppanteI; Annalisa De PaolisI; Raffaella PecciIII; Rossella BediniIII
IDipartimento di Ingegneria Meccanica e Aerospaziale, Sapienza Università di Roma, Rome, Italy
IIOrtopedia e Traumatologia, Policlinico Universitario Campus Bio-Medico, Rome, Italy
IIIDipartimento di Tecnologie e Salute, Istituto Superiore di Sanità, Rome, Italy
ABSTRACT
INTRODUCTION: Micro-CT analysis is a powerful technique for a non-invasive evaluation of the morphometric parameters of trabecular bone samples. This elaboration requires a previous binarization of the images. A problem which arises from the binarization process is the partial volume artifact. Voxels at the external surface of the sample can contain both bone and air so thresholding operates an incorrect estimation of volume occupied by the two materials.
AIM: The aim of this study is the extraction of bone volumetric information directly from the image histograms, by fitting them with a suitable set of functions.
METHODS: Nineteen trabecular bone samples were extracted from femoral heads of eight patients subject to a hip arthroplasty surgery. Trabecular bone samples were acquired using micro-CT Scanner. Hystograms of the acquired images were computed and fitted by Gaussian-like functions accounting for: a) gray levels produced by the bone x-ray absorption, b) the portions of the image occupied by air and c) voxels that contain a mixture of bone and air. This latter contribution can be considered such as an estimation of the partial volume effect.
RESULTS: The comparison of the proposed technique to the bone volumes measured by a reference instrument such as by a helium pycnometer show the method as a good way for an accurate bone volume calculation of trabecular bone samples.
Key words: microtomography, osteoarthritis, osteoporosis, bone morphometry, image processing
INTRODUCTION
Morphometric characterization of trabecular bone samples is an open research topic since it can be useful for the study and early diagnosis of pathologies such as osteoarthritis or osteoporosis [1-5]. This characterization allows to analyze the bone structure and it is also helpful for the mechanical characterization of bone since the mechanical properties are usually related to morphometric parameters.
Histology is the standard procedure for morphometric analysis of bone. Although it is a good technique, it traditionally requires a destructive preparation process (i.e. sectioning, staining). Therefore, in the last years many studies worked on the bone characterization by an alternative technique: micro-Computed Tomography (micro-CT) [6-10]. This technique allows to obtain an accurate visualization and quantification of cancellous bone microstructure and furthermore, differently from the traditional histology, it does not require the destructive preparation of the sample.
One of the main problem of bone analysis by micro- CT is the processing of the scans for sample reconstruction. This usually requires a process named binarization which consists on the definition of a threshold value of grey-level, necessary to distinguish the bone from the background.
The choice of this value is a crucial task since a standard method doesn't exist, so it has been the subject of several studies and different solutions [11]. Ding et al. [12] calculates threshold value by an adaptive method which consists on the segmentation of the grey levels data at different thresholds. The correct threshold is then chosen considering where the volume fraction changes the least. Other parameters can be also considered to define the threshold value. For example, observing that in a certain range of grey levels, density connectivity has a constant value and structure model index is linear, the mean value of this range can be chosen as threshold value [13]. In a previous study [14], a different method was implemented in order to calculate threshold value of bone samples. This method is based on the histogram analysis and requires the splitting of the histogram in three different populations: pixels marked as bone, pixels marked as air and a third population of pixels that can not be statistically referred to bone or air. A statistical analysis was then used in order to define a threshold value.
The choice of the correct threshold value is not the only goal of this study, the resolution spatial sampling of bone elaboration of micro-CT system presents another problem, related to the fact that bone is an inhomogeneous material. From the binarization process, each voxel is associated to bone or air, not considering that each voxel of the acquired image can be composed by both of them. This effect is called "partial volume effect" (PVE) and it affects especially voxels at the edges of the analyzed sample [15-19].
The aim of this study is the application of a technique [20] for the analysis of the image histograms of the sample which allows the classification of the voxels into a) classes relative to the different materials that are present into the imaged volume and b) classes relative to the voxels that are composed by a "mixture" of two different materials. Without any further data processing, the actual volume that should be attributed to a single material (bone in our case) within these latter classes of voxels would be unknown and this circumstance causes the so called PVE. Moreover, this technique allows to overcome the problems related to the binarization process, which can introduce significant errors when the actual bone volume of a sample has to be measured. In fact, depending on the chosen threshold, many morphometric parameters, such as the bone volume, can be under or overestimated.
The histogram of a homogeneous material is composed of a peak centered in the grey level correspondent to the material so the volume of the sample could be calculated as the product between voxel volume and number of bone voxels.
In actual conditions, bone inhomogeneity implies that bone voxels can be represented by more than one grey level. Because of this, the grey level histogram of a bone image is composed of two ideal pulses (considering air as a second material) spread over the entire range of grey levels.
In a previous work [14], we carried out a detailed segmentation of the histograms of micro-CT trabecular bone images. From this analysis emerges that the histograms can be decomposed into four different classes, estimated by Gaussian functions. These components correspond to voxels that represent 1) bone, 2) intertrabecular voids, 3) the void space outside the sample and 4) the voxels that lies across the bone/void boundary.
As referred above, the present work starts from a segmentation method proposed by Laidlaw [20] for MRI images.
According with this approach, each histogram can be split into two base functions corresponding to the bone and void portions of the sample and a third function representing voxel that contain a mixture of bone and air. This third contribution can be considered as an estimation of the PVE. Under these hypotheses and comparing the obtained results with the bone volume of each sample measured by a helium pycnometer, we found that the bone volume can be computed as the sum of voxels in the bone distribution plus the 50% of voxels of the third distribution.
MATERIALS AND METHODS
Data acquisition
Nineteen trabecular bone samples were extracted from femoral heads of eight patients subject to a hip arthroplasty surgery. The samples were collected as waste from the operating room and thus destined to destruction. After being catalogued, samples were preserved in a freezer at 10 ºC for about a month. From the middle of the femur capita of each sample was obtained a slice of about 10 mm of thickness corresponding to the frontal plane and then these slices were kept 10 hours in a freezer.
Subsequently, slices were subjected to three complete cycles of dehydration with aqueous solutions having crescent percentage of ethanol, severally 70%, 90% and 99.9%, in order to defatted them. Between the dehydratation cycles, samples were kept in a refrigerator for about 10 hours. Slices were then cut with a diamond belt saw (EXTEC Labcut 1010, Enfield CT) in order to obtain cubic samples of about 6 mm. These samples were further dehydrated and defatted by other three cycles with ethanol as previously described [21, 22].
Trabecular bone samples were acquired using a Skyscan 1072® micro-CT Scanner [23]. The acquisition parameters were set at 100 kV and 98 µA, with 1 mm of aluminum filter and with an isotropic voxel size in the range of 11.24 µm-14.66 µm. An angular step of 0.45 degree was used and 400 projections were acquired over an angular range of 180 degrees. Using Cone Beam Reconstruction, software based on the Feldkamp algorithm [24], from these projection images in TIFF format, transversal sections of the sample in BMP format can be obtained (Figure 1). The distance between two consecutive sections is the double of the resolution chosen for the acquisition.
After the reconstruction, data were processed by a software of Skyscan®, named Ct-Analyzer®. This software allows to obtain histomorphometric parameters of the samples, such as those described in [25].
As described above, because of bone inhomogeneity, noise and the finite resolution of the CT apparatus, each peak of the histogram was spread into a Gaussian-like function. The first peak is centered on the mean gray level produced by the bone x-ray absorption, the second one represents the portions of the image occupied by air. In the middle of these two functions, a third distribution, representing voxels that contain a mixture of bone and air, can be detected. This contribution can be considered such as an estimation of the partial volume effect.
In order to verify the accuracy of the implemented method, the volume of the bone samples were also measured by a helium pycnometer (AccuPyc II 1340, Micromeritics-GA, USA) [26].
Data processing
From the images of the transversal sections of the sample, the histogram of each sample was computed.
Then [20], the histogram can be fit with two Gaussian basis functions:
where g represents the histogram grey levels, c and s represent mean and standard deviation, respectively. Finally gi, ci, siare scalar component of g, c, s. The histogram component relative to the mixture along boundaries between two materials (bone and air, in our case) is calculated with this third basis function:
An algorithm was implemented in order to fit the histogram obtained by the acquired images, starting from the equations of the three basis functions fbone, fair, fmixture. The fitting process was based on the Levenberg-Marquardt algorithm [25], which is an iterative method to solve non-linear least squares problems.
In order to compute the bone volume to assign to the fmixturecomponent of the histogram, we made the realistic assumption: the gray level g of the voxels of the mixture is due to the portion of bone within the voxel itself, represented by the random variable bv, with 0 < bvi <1. This is clearly due to the circumstance that only the portion bv contributes to the attenuation of x-rays.
We can express the bone volume of the mixture BVmixtureas the integral of the product between bv and its probability density function, multiplied by the overall volume of the voxels pertaining to the mixture Vmixture, considering also that, according to [20] the mean value of bv is 0.5 (i.e. one half of the voxel volume).
This result clearly does not depend on the local bone density. Subsequently, we can estimate Vmixtureusing fmixture (g, c1, c2, s):
where Vv is the volume of a voxel. In the end we obtain:
We can conclude that the bone volume could be assessed as the sum of the volume of all the voxels pertaining to fbone, plus the volume of 50% of the voxels that belong to fmixture.
Volumetric data were then compared with the volume computed with the reference method. The helium pycnometer, in our measuring conditions, allowed to determine the sample volume with a reproducibility typically of ± 3 mm3 of the full-scale cell chamber volume, i.e. 10 cm3.
RESULTS
According to equations (1) and (2), the implemented algorithm allows to fit each histogram with two Gaussian distributions and a distribution relative to the mixture. A typical histogram and the distributions for one of the analyzed samples are reported in Figure 2a.
As described above, volume of the samples is calculated as the sum of voxels in the bone distribution plus the 50% of voxels of the third distribution. Volumes computed by this method are reported in Table 1, together with those measured by the pycnometer (mean and standard deviation).
Comparing these values, an excellent correlation can be established, as reported in the graph of Figure 2b.
Angular coefficient of the linear fit, equal to 0.9854, shows a slight underestimation of computed vs measured volume. The intercept of the linear fit, equal to about 1.3 mm3, lies within the accuracy of the pycnometer, that is, ± 3 mm3 and ± 0.03% of the reading. These values have been compared with the percentage error calculated as the difference between values computed with the proposed method and the reference values measured by the pycnometer.
Finally, the viability of the implemented method has been verified comparing the calculated volumes with those resulting from a previous study, based on the implementation of a previously proposed method for the image binarization [14].
An example of this comparison is reported in the graph in Figure 2c in which is evident that even small differences of binarization threshold will lead to great variations of the computed volume. This trend is depicted with the black continuous line.
Considering the value obtained by pycnometer as reference, it can be observed that the volume computed with this new method is more accurate than the previous method. Both methods focus on the histogram analysis, however the first method leads to define a threshold for the image processing, this new method avoid the problem of the threshold definition because it allows to extract bone volumetric information directly from the histogram.
CONCLUSION
Bone characterization is a crucial research issue because of its importance in the prediction of diseases such as osteoarthritis or osteoporosis.
Because of the importance of this topic, in recent years many studies worked on the bone morphometric characterization by microtomographic technique. It is a good investigation tool since, instead of traditional technique such as histology, it allows an accurate but non invasive analysis. The calculation of the parameters is based on the analysis of the images which requires the definition of a proper grey level threshold value. The definition of the threshold is a crucial point because, despite this topic has been the subject of many studies, a reliable procedure to define it is not available yet.
The aim of this study is to test a new technique for the bone volume calculation in order to overcome the threshold problem in micro-CT images, and more generally in all the images that are affected by PVE.
Nineteen samples of bone cube, extracted from femoral heads of eight patients subject to a hip arthroplasty surgery, were acquired by a micro-CT apparatus. The implemented algorithm extrapolates bone volumetric information by the analysis of the grey level histogram of the reconstructed images. Volume values of the samples were also measured with a helium pycnometer, taken as a reference.
A good correlation between the two datasets were demonstrated so the histogram processing can be defined a viable method for an accurate and non destructive bone volume calculation.
Acknowledgements
The authors wish to thank P. Scarlato and S. Mollo from Italian National Institute of Geophysics and Volcanology (INGV) for providing the helium pycnometer and for their assistance during the measurements.
Conflict of interest statement
There are no potential conflicts of interest or any financial or personal relationships with other people or organizations that could inappropriately bias conduct and findings of this study.
REFERENCES
1. Legrand E, Chappard D, Basle MF, Audran M. Evaluation of trabecular microarchitecture. Prospects for predicting the risk of osteoporosis and fracture. Rev Rhum [Engl Ed] 1999;66:543-7.
2. Dougherty G. Quantitative CT in the measurement of bone quantity and bone quality for assessing osteoporosis. Med Eng Phys1996;7:557-68. DOI: 10.1016/1350- 4533(96)00011-2
3. Ding M, Odgaard A, Hvid I. Changes in the three dimensional microstructure of human tibial cancellous bone in early osteoarthritis. J Bone Jt Surg Br 2003;85:906-12.
4. Fazzalari NL, Parkinson IH. Femoral trabecular bone of osteoarthritic and normal subjects in an age and sex matched group. Osteoarthritis and Cartilage 1998;6:377- 82. DOI: 10.1053/joca.1998.0141
5. Parfitt AM, Mathews CHE, Villanueva AR, Kleerekoper M, Frame B, Rao DS. Relationships between Surface, Volume, and thickness of iliac trabecular bone in aging and in osteoporosis. Implications for the microanatomic and cellular mechanisms of bone loss. J Clin Invest 1983;72(4):1396-409. DOI: 10.1172/JCI111096
6. Kuhn JL, Goldstein SA, Feldkamp LA, Goulet RW, Jesion G. Evaluation of a microcomputed tomography system to study trabecular bone structure. J Orthop Res 1990;8(6):833-42. DOI: 10.1002/jor.1100080608
7. Odgaard A. Three-dimensional methods for quantification of cancellous bone architecture. Bone 1997;20(4):315-28. DOI: 10.1016/S8756-3282(97)00007-0
8. Jones CA, Arns CH, Sheppard AP, Hutmacher DW, Milthorpe BK, Knackstedt MA. Assessment of bone ingrowth into porous biomaterials using MICRO-CT. Biomaterials 2007;28(15):2491-504. DOI: 10.1016/j.biomaterials.2007.01.046
9. Muller R, Ruegsegger P. Micro-tomographic imaging for the non-destructive evaluation of trabecular bone architecture. Stud Health Technol Inform 1997;40:61-79.
10. Hildebrand0 T, Laib A, Muller R, Dequeker J, Ruegsegger P. Direct three dimensional morphometric analysis of human cancellous bone: microstructural data from spine, femur, iliac crest, and calcaneus. J Bone Miner Res 1999;14:1167-74.
11. Tan Y, Kiekens K, Kruth J, Voet A, Dewulf W. Material dependent thresholding for dimensional x-ray computed tomography. In: Tan Y (Ed.). DGZFP-Proceedings BB 128-CD. International Symposium on Digital Industrial Radiology and Computed Tomography, Berlin, Germany, 20-22 June 2011.
12. Ding M, Odgaard A, Hvid I. Accuracy of cancellous bone volume fraction measured by micro-CT scanning. J Biomechanics 1999;32(3):323-6. DOI: 10.1016/S0021- 9290(98)00176-6
13. Kim CH, Zhang H, Mikhail G, von Stechow D, Muller R, Kim HS, Guo XE. Effects of thresholding techniques on micro-CT-based finite elemet model of trabecular bone. J Biomech Eng 2007;129(4):481-6. DOI: 10.1115/1.2746368
14. Bedini R, Marinozzi F, Pecci R, Angeloni L, Zuppante F, Bini F, Marinozzi A. Analisi microtomografica del tessuto osseo trabecolare: influenza della soglia di binarizzazione sul calcolo dei parametri istomorfometrici. Roma: Istituto Superiore di Sanità; 2010 (Rapporti ISTISAN, 10/15).
15. Link TM, Majumdar S, Grampp S, Guglielmi G, van Kuijk C, Imhof H, Glueer C, Adams JE. Imaging of trabecular bone structure in osteoporosis. Eur Radiol 1999;9:1781-8. DOI: 10.1007/s003300050922
16. Muller R, Koller B, Hildebrand T, Laib A, Gianolini S, Ruegsegger P. Resolution dependency of microstructural properties of cancellous bone based on threedimensional mu-tomography. Technol Health Care 1996;4:113-9.
17. Peyrin F, Salome M, Cloetens P, Laval-Jeantet AM, Ritman E, Ruegsegger P. Micro-CT examinations of trabecular bone samples at different resolutions: 14, 7 and 2 micron level. Technol Health Care 1998;6:391-401
18. Kothari M, Keaveny TM, Lin JC, Newitt DC, Genant HK, Majumdar S. Impact of spatial resolution on the prediction of trabecular architecture parameters. Bone 1998;22:437-43. DOI: 10.1016/S8756-3282(98)00031-3
19. Tabor Z. Analysis of the influence of image resolution on the discriminating power of trabecular bone architectural parameters. Bone 2004;34:170-9. DOI: 10.1016/j. bone.2003.10.002
20. Laidlaw DH, Fleisher KW, Barr AH. Partial volume Bayesian classification of material mixtures in MR volume data using voxel histograms. IEEE Transaction Medical Imaging 1998;17(1):74-86. DOI: 10.1109/42.668696
21. Marinozzi F, Marinozzi A, Bini F, Zuppante F, Pecci R, Bedini R. Variability of morphometric parameters of human trabecular tissue from coxo-arthritis and osteoporotic samples. Ann Ist Super Sanità 2012;48:19-25. DOI: 10.4415/ANN_12_01_04
22. Marinozzi F, Zuppante F, Bini F, Marinozzi A, Pecci R, Bedini R. Bone volume measurement of human femur head specimens by modeling the histograms of micro-CT images. In: Di Gianberardino et al. (Eds). Computational modelling of objects represented in images: Fundamentals, methods and applications III. Leiden, The Netherlands: CRC Press/Balkema; 2012.
23. Skyscan. 1072 X-ray microscope instruction manual. Aartselaar (Belgium): Skyscan Ed.; 2005.
24. Feldkamp LA, Davis LC, Kress JW. Practical cone-beam algorithm. JOSA 1984;1(6):612-9. DOI: 10.1364/ JOSAA.1.000612
25. Parfitt A, Drezner K, Glorieux H, Kanis A, Malluche H, Meunier PJ, Ott SM, Recker RR. Bone histomorphometry: Standardization of nomenclature, symbols and units. JBMR 1987;2(6):595-610. DOI: 10.1002/ jbmr.5650020617
26. Tamari S, Aguilar-Chávez A. Optimum design of the variable-volume gas pycnometer for determining the volume of solid particles. Meas Sci Technol 2004;15:1146-52. DOI: 10.1088/0957-0233/15/6/015
Address for correspondence:
Franco Marinozzi
Dipartimento di Ingegneria Meccanica e Aerospaziale, Sapienza Università di Roma
Via Eudossiana 18
00184 Rome, Italy
E-mail: franco.marinozzi@uniroma1.it
Received on 21 December 2012.
Accepted on 16 July 2013.