An Independent Active Contours Segmentation framework for bone micro-CT images

Computers in Biology and Medicine 87 (2017) 358–370

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/compbiomed

An Independent Active Contours Segmentation framework for bone micro-CT images Vasileios Ch. Korfiatis a, Simone Tassani b, George K. Matsopoulos a, * a b

School of Electrical and Computer Engineering, National Technical University of Athens, Greece Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain

A R T I C L E I N F O

A B S T R A C T

Keywords: Micro-CT Trabecular bone Image segmentation Active Contours ROI extraction

Micro-CT is an imaging technique for small tissues and objects that is gaining increased popularity especially as a pre-clinical application. Nevertheless, there is no well-established micro-CT segmentation method, while typical procedures lack sophistication and frequently require a degree of manual intervention, leading to errors and subjective results. To address these issues, a novel segmentation framework, called Independent Active Contours Segmentation (IACS), is proposed in this paper. The proposed IACS is based on two autonomous modules, namely automatic ROI extraction and IAC Evolution, which segments the ROI image using multiple Active Contours that evolve simultaneously and independently of one another. The proposed method is applied on a Phantom dataset and on real datasets. It is tested against several established segmentation methods that include Adaptive Thresholding, Otsu Thresholding, Region Growing, Chan-Vese (CV) AC, Geodesic AC and Automatic Local RatioCV AC, both qualitatively and quantitatively. The results prove its superior performance in terms of object identification capability, accuracy and robustness, under normal circumstances and under four types of artificially introduced noise. These enhancements can lead to more reliable analysis, better diagnostic procedures and treatment evaluation of several bone-related pathologies, and to the facilitation and further advancement of bone research.

1. Introduction Microtomography (micro-CT) using X-rays is a scanning method for the high-resolution imaging of living tissues (in vivo or ex vivo) and objects of few micrometers, which has gained increased popularity especially as a pre-clinical application [1,2]. It permits the threedimensional (3D), non-destructive investigation of a specimen with relatively low cost and scanning efficiency, without any special specimen preparation [3]. One of the major applications of micro-CT scanning in living organisms is bone analysis for the investigation of various bone pathologies, such as osteoarthritis [4] and osteoporosis [5], bone functions, such as bone remodeling [6], biomechanical properties [7] and response to treatment [8,9]. Rodents are the most common samples, with typical in-vivo scanning spatial resolution that of 50 μm [10,11]. The study of trabecular bone requires higher resolution due to its thickness, and is often performed ex vivo using excised bones. This also applies to human bones [3,12]. Bone analysis is performed by extracting several morphometrical parameters on the binary image obtained by segmenting the original micro-CT volume [13,14]. The performance of the

* Corresponding author. 9, Iroon Polytechniou str, Zografos, 15780, Athens, Greece. E-mail address: [email protected] (G.K. Matsopoulos). http://dx.doi.org/10.1016/j.compbiomed.2017.06.016 Received 10 November 2016; Received in revised form 16 June 2017; Accepted 16 June 2017 0010-4825/© 2017 Elsevier Ltd. All rights reserved.

segmentation step is critical due to the high sensitivity of the morphometrical parameters to errors/deviations [15]. Currently, there is no accepted and well-established method for micro-CT segmentation [1,16,17]. The most widely used algorithmic segmentation methods are based on intensity histogram thresholding, either by a manual selection of a global threshold [1,3,12,18,19] or by semi-manual selection using an initial human estimation followed by algorithmic refinement [20]. Although the above approaches are easy to implement, they have showed to introduce significant errors when variable bone structures are involved [21] while their need for manual intervention makes them time-consuming for the experts and subjective. Other intensity threshold-based segmentation techniques have been proposed, including Otsu Thresholding [22], Hysteresis Thresholding [23], Dual-Thresholding [24] and Adaptive Thresholding [25]. However, these methods remain sensitive to noise, provide medium accuracy (especially on high resolution images), still require some form of manual intervention and may miss some important parts of the foreground, such as smaller but relevant trabeculae. In order to overcome these limitations, Region Growing (RG) [26], Ray Casting Algorithm (RCA) [27] and

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

Active Contours [28,29] techniques have been recently proposed as alternatives. Active Contours (AC) [30] are variational segmentation methods that are based on the theory of curve evolution. Implicit AC models substitute the explicit curve representation with a Level Set (LS) function [31], and have been used successfully for the segmentation of synthetic and real images [32,33]. Implicit AC are categorized into edge-based [30,34], or region-based [35,36], and global [34–36] or local [32,33,37]. Chan-Vese (CV) [35] is an implicit global AC algorithm that aims at minimizing the intra-region variance, by extending the Mumford-Shah energy functional [38]. The CV AC enables multiple objects' segmentation, providing higher accuracy and object identification ability. However, it is computationally heavy, sensitive to initialization and needs manual determination of its optimal parameters' values. These drawbacks lead to medium robustness and slower execution times. Modifications based on local information may provide better results, but increase algorithmic complexity, execution times per iteration and number of manually defined parameters. In this paper, a segmentation framework for two dimensional (2D) trabecular bone micro-CT images, called Independent Active Contours Segmentation (IACS) is presented, focusing on the identification of single trabeculae, which can appear much smaller than typical trabeculae and thus easy to miss during typical segmentation schemes. It is comprised of two main modules: Automatic Region of Interest (ROI) Extraction and Independent AC (IAC) Evolution. The automatic ROI extraction is a combination of established unimodal thresholding concepts with a novel heuristic thresholding technique for noise reduction of binary images. The Independent AC Segmentation algorithm uses a modified version of the CV AC algorithm, where an independent AC instance is assigned to each extracted ROI and evolves strictly within it. To the author's knowledge, IACS introduces five novelties, namely IAC Evolution, which enables the use of multiple independent ACs in the same image, the Automatic ROI Extraction module, which can also be used for AC initialization, the Object Maintaining Stepwise Thresholding technique, which serves as a binary image denoising process, the AutoCoefficient modification for the automatic calculation of the CV AC force coefficients and finally, the successful application of all the above independent features in a newly unified framework for the automatic microCT data segmentation. The proposed scheme is evaluated on a micro-CT Phantom dataset [2] and on real trabecular bone micro-CT images. For the Phantom dataset, three-dimensional (3D) reconstructions are used for the qualitative evaluation, and the original dataset along with noise-distorted versions (using four noise types) are used for the quantitative evaluation. For the real images, qualitative evaluation is performed using 3D reconstructions and 2D case studies. The proposed methodology is compared against several well-established segmentation techniques, such as Otsu, Adaptive Thresholding, and popular AC algorithms, providing a highly accurate, robust, low-noise and automatic framework for trabecular bone segmentation. This paper is structured as follows: in Section 2 the proposed IACS framework is introduced, along with a brief description of the rest of the algorithms under comparison. In Section 3, the Phantom dataset and the methodology for creating its optimized segmentation equivalent are presented, followed by the qualitative and quantitative results from the Phantom dataset and real data, illustrating the abilities of the proposed IACS framework. Finally, on Section 4, further discussion is presented on the capabilities of the proposed segmentation framework and on specific selected relevant concepts of the methodology.

Fig. 1. Block Diagram of the proposed IACS framework.

foreground. This coarse pre-segmentation serves as input for IAC Evolution module, which applies one independent instance of the modified CV AC segmentation algorithm in each ROI and returns the final segmentation result. 2.1. Automatic ROI extraction The Automatic ROI Extraction module is comprised of four steps: Background Elimination, Object Maintaining Stepwise Thresholding, Image Patching and ROI Extraction. 2.1.1. Background Elimination step The Background Elimination (BE) step is based on the Triangular histogram thresholding method, which has been successfully used in the past for unimodal histogram segmentation [39–41]. As it can be seen in Fig. 2(b), the intensity histogram of a typical micro-CT trabecular bone image, shown in Fig. 2(a), can be safely considered unimodal. Specifically, the background corresponds to the central part of the distribution, which is located on the higher histogram values and the foreground corresponds to the long tail of the distribution, which extends to the left to lower histogram values. It should be noted here that in the given datasets higher attenuation objects have lower intensity, thus the actual objects appear darker (lower part of the histogram) and the background whiter, which is the opposite of typical CT images. Triangular histogram thresholding draws a line from the peak of the distribution to the last non-zero bin and then calculates the vertical Euclidean distance of all the bins from this line. The bin with the largest distance from the line is chosen as the threshold level. The image then is transformed according to:

argvalðIðx; yÞ  ThresholdÞ→Iðx; yÞ ¼ Threshold

(1)

This process has one user-defined parameter, the number of iterations, which depends on the number of backgrounds present in the image. This parameter definition is justified by the fact that micro-CT images contain multiple and generally homogenous backgrounds. Each background has a large central histogram value and a small standard deviation. Moreover, the central values of the backgrounds' distributions are located very close to each other. Consequently, the individual background distributions overlap to a great degree, formulating the visible peak of the actual histogram distribution (Fig. 2(b)). By using the Triangular histogram thresholding these individual background components are removed, one at a time, and the foreground details emerge more clearly in the histogram. Current micro-CT images contain three backgrounds, which correspond (from the outside to the inside) to the air surrounding the sample-holder, the plastic of the sample-holder itself and the water inside the holder where the trabecular bone specimen is submerged. Consequently the process is repeated three times. The final result is illustrated in Fig. 2(c).

2. Methodology The proposed IACS framework, illustrated in Fig. 1, is comprised of two modules: Automatic ROI Extraction and IAC Evolution. The Automatic ROI Extraction module is a pre-processing step aiming at removing the majority of image's background, while maintaining all the 359

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

Fig. 2. Illustration of the Automatic ROI Extraction steps. (a) Original Image, (b) Histogram of (a), (c) Result of Background Elimination step, (d) Result of Object Maintaining Stepwise Thresholding, (e) Patched image, (f) Extracted ROI.

from the contour of the object to its center (where it will have the lowest values in the current representation) and that the standard deviation of the actual foreground object's intensity is not small. This is not true for the micro-CT backgrounds, which do not present a constant increase or decrease as we move in any direction, due to their specific textural properties, and present a small intensity standard deviation. Consequently, as the threshold decreases for a sufficient number of repetitions foreground objects will shrink but will be maintained as entities, but the remaining background information or artifact noise will not, with very high probability. As a result, the number of objects will generally decline as the number of iterations increases, until it reaches a stable solution. An exception of this behavior is when a large object may give rise to multiple objects, due its separation. On that case the number of objects will

2.1.2. Object Maintaining Stepwise Thresholding step After BE step, a novel methodology called Object Maintaining Stepwise Thresholding (OMST) is applied to the transformed image. This is an iterative denoising process that removes large remaining noise and background artifacts. The method is initialized using the final threshold of the BE step. On every iteration, it segments the image using the current threshold, finds the number of objects, decreases the threshold value by one and checks if the number of objects has changed. The process ends when there is no change in the number of objects for two consequent iterations. The result is a binary image that contains a coarse segmentation of the foreground, seen in Fig. 2(d). The concept behind this step is based on the observations that on actual foreground objects the intensity steadily decreases as we move 360

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

corresponding ACs level set functions. The IAC evolution is as follows:

increase but only for the current iteration; then it should stay the same for the rest of the algorithm's execution as the aforementioned behavior takes place.

1. Calculate the c1;j and c2; j intensity means inside and outside φj , for all the pixels inside ROI j. 2. Calculate the evolution forces Forcej ðx; yÞ for each ROI j using Eq. (3), according to the chosen AC data term (CV in our case). 3. Set Forceðx; yÞ ¼ 0 for the rest of the image (which does not belong to a ROI). 4. Evolve φj according to Step 2. 5. Regularize φj using Gaussian filtering. 6. Repeat Steps 1–5 until the maximum number of stable iterations endingIter for φj has been reached.

2.1.3. Image Patching and ROI extraction steps The result of OMST step is patched using fixed size patches, for the Image Patching step (Fig. 2 (e)), and the patches that contain foreground pixels are identified and connected into ROIs using 4-connectivity, for the ROI extraction step (Fig. 2 (f)). The patch size is a user-defined parameter, tuned in order to provide balance between ROIs' detail level and maintained background information. This parameter's value depends on the image size and an illustration of its influence on the ROI extraction step can be found in Fig. 3, where the results using 8  8 (Fig. 3 (a)), 16  16 (Fig. 3 (b)), 32  32 (Fig. 3 (c)) and 64  64 (Fig. 3 (d)) patch sizes are presented, with different colors highlighting the identified ROIs on each case. It is evident that smaller patch sizes translate to smaller ROIs with less background information, with the cost of more computational time.

The above process takes place simultaneously for all φj . The intensity means, evolution forces and current number of stable iterations (stopping criterion) are independent for each φj . As a result, the whole AC evolution is tailored to the local characteristics of each ROI. The forces outside the boundaries of each ROI are set to zero on every iteration in order to remove the possibility that the AC will evolve outside of them. Regarding the CV AC data term modification, an automatic scheme is proposed according to the following function:

2.2. Independent Active Contours evolution In the IAC Evolution module, each ROI is assigned a dedicated AC instance, which evolves independently of the others according to the ROI's local attributes. The AC used in this paper is a novel modified version of the Chan-Vese AC algorithm [35] that automatically assigns values to the force coefficients (λ1 and λ2) according to the intensity difference between the inside and the outside of the AC. As a result, it removes two user-defined parameters, namely the force coefficients. Let j ¼ 1; …; N the number of ROIs in the grey-level image and φj the

ðmaxIntensityj c2;j Þ=jc1;j c2;j j AutoCoefficient ¼ e

(2)

and the proposed CV AC data term is now:


2
2 Forcej ðx; yÞ ¼ AutoCoefficient⋅ Iðx; yÞ  c1;j  Iðx; yÞ  c2;j

(3)

Fig. 3. Influence of patch size on a real micro-CT image using a (a) 8  8, (b) 16  16, (c) 32  32 and (d) 64  64 patch size. The final ROIs are highlighted using various colors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 361

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

2. Adaptive Thresholding (AT): Adaptive Thresholding (AT) is a local thresholding method that returns different thresholds for each image neighborhood. It has three main parameters, neighborhood size, mean intensity calculation method (actual mean or weighted mean) and subtraction constant from mean [25]. 3. Region Growing (RG): Region Growing (RG) is an iterative region based segmentation method that adds neighboring pixels of current foreground points according to a similarity criterion, usually intensity value. Besides the initial seeds mask, it has one parameter, the acceptable max-min intensity difference [26]. 4. Original Chan-Vese Active Contour (CV AC) [35]. 5. Automatic Local Ratio Chan-Vese Active Contour (ALR-CV) [37]. 6. Geodesic Active Contour (GAC) [34].

where c1;j is the mean intensity inside φj , c2;j is the mean intensity outside φj , Iðx; yÞ is the intensity of pixel in (x,y) position, and maxIntensityj is the maximum intensity of the ROI. This modification automatically   increases the outside forces when the c1;j  c2;j  (foreground-background) difference is small, allowing for tighter contours around fainter dark objects. It also has the ability to change its value on every iteration, adapting to the current requirements of the algorithm. Following typical level-set evolution schemes, the Forcej ðx; yÞ term is normalized and multiplied by the time step term dt ¼ 0:5, and then used as a magnitude for the zero level-set curve evolution in the normal direction [30,31,35]. The normalization is done separately for each ROI j by dividing each Forcej ðx; yÞ with its absolute maximum value. Eq. (3) has negative values inside of the contours and positive outside of the contours. The proposed AutoCoefficient modification, compared to the original algorithm, enables the further adaptation of the proposed scheme on the various parts of the image, resulting in higher overall accuracy, while being a completely automatic procedure. The effect of the proposed CV AC AutoCoefficient modification on the segmentation result is illustrated in Fig. 4. Specifically, the image contains four objects, ranging from very clearly defined (green) to very faint (purple), while the segmentation result of the CV AC without the modification, using λ1 ¼ λ2, is shown in Fig. 4 (a) and the corresponding result after the use of the proposed modification is shown in Fig. 4 (b). Three out of the total four objects (green, red and light blue) are segmented similarly in both images. However, for the faintest fourth object (purple), the original version returns a segmentation area similar to that of the other three objects (Fig. 4 (a)), which can be characterized as an oversegmentation, whereas the modified version returns a clearly smaller segmentation area (Fig. 4 (b)), which constitutes a better approximation. For the initialization mask, the darkest 1% of pixels of each ROI are used as seeds, since the foreground is darker than the background. For the regularization, Gaussian filtering replaced the typical regularization term resulting in a faster process [42,43], a modification that has been used in many works [32,33,37,44]. The regularization parameter r is set to 0.25 for all the AC schemes, including the proposed one, and the Gaussian filter kernel has size 5  5.

All the Active Contours (CV, ALR-CV and GAC) have a contour regularization parameter r and an initialization mask, but CV and GAC have two more parameters, the values of force coefficients λ1 and λ2 . For the AT, the calculation method used was actual mean, neighborhood size was set to 17 and local threshold mean difference to 31. RG was initialized with seeds inside every object (manual annotation) and max-min intensity difference was set to 20. All the ACs have a contour regularization parameter r ¼ 0:25, λ1 ¼ 1:6 and λ2 ¼ 1, are initialized using the typical scheme with a circle of Radius ¼ 100 px in the center of the image (where the foreground exists), and have maxIterations ¼ 10000 and endingIter ¼ 3. The same values for maxIterations and endingIter are also used in the proposed IACS scheme. The above configurations were chosen in order to provide tighter contours and smaller areas. 3. Results 3.1. Data acquisition The micro-CT Phantom is comprised of 991 gray-level micro-CT slices with intensity values in the range [0, 255], a size of 1024  1024 pixels and an isotropic voxel size of 19.47 μm (model Skyscan 1072, Bruker microCT, Kontich, Belgium) [2]. The dataset contains 15 objects of various sizes and shapes, whose thickness ranges from 20 μm to 1000 μm (1 μm ¼ 106 m). It has been found to provide a very good approximation of the real-case scanning environment and a guaranteed thickness value by the manufacturer, which enables the analysis of segmentation errors due to micro-CT measurements. Finally, the images are inverted, meaning that higher attenuation objects have lower intensity.

2.3. Comparison with existing micro-CT segmentation algorithms Based on the current state of the literature on bone micro-CT segmentation and AC, six algorithms were implemented and tested in order to compare their results with those of the proposed IACS. These are:

3.2. Micro-CT Optimized Phantom dataset creation

1. Otsu Thresholding: Otsu is a thresholding algorithm that automatically determines the optimal threshold value of the intensity histogram. As a result, it does not have any parameters [22].

In previous works using the aforementioned Phantom dataset, the evaluation of the segmentation results was made using only the morphometrical parameter of object's volume thickness. In this work, an

Fig. 4. Influence of the proposed AutoCoefficient modification of the Chan-Vese algorithm on the segmentation result of a Phantom micro-CT image. (a) Original Chan-Vese Algorithm, (b) Modified Chan-Vese Algorithm with AutoCoefficient. 362

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

optimized segmentation dataset for the Phantom micro-CT is introduced for the first time, called hereafter micro-CT Optimized Phantom (OP). This novelty enables the use of established image processing evaluation metrics directly on the segmentation results' images, allowing for a more technical evaluation approach. The micro-CT OP calculation was based on object's thickness and dimensions, pixel size information (pixel size ¼ 19.47 μm) and the use of “CTAnalyzer” software, Bruker [45]. . Specifically, each object was isolated by cropping the corresponding slices and its sub-dataset was used as input to CTAnalyzer. Its volume thickness was calculated for all thresholds from 0 to 255, keeping the threshold value that provided the smallest deviation from the object's nominal thickness. Finally, the cropped segmented patch was placed on the OP mask, in its original coordinates. The same procedure was followed exactly in the seminal paper [2], although the authors did not use cropped sub-images for each object but proposed a global threshold that minimized the total volume thickness error. The two smallest wires (20 μm and 50 μm) and the smallest foil (20 μm) could not be standardized using the above procedure due to their small size. The 50 μm wire and the 20 μm foil, were manually placed in the OP using the reported pixel size ¼ 19.47 μm and the objects' dimensions and coordinates. The 20 μm wire was practically invisible even by human experts, as a result it was omitted from the micro-CT OP. The optimal threshold values along with the corresponding measurements are reported in Table 1. It is clear that micro-CT OP provides a much better object approximation of the Phantom than that of the FT-141, while enabling the use of technical evaluation measures.

Index of Balanced Harmonic F Score ¼ ð1 þ a⋅ðSensitivity  SpecificityÞÞ⋅Harmonic F Score

Index of Balanced GMean ¼ ð1 þ a⋅ðSensitivity  SpecificityÞÞ⋅GMean (7) where a ¼ 1, is the dominance parameter. Sensitivity (or Recall) equals true positive rate, which indicates how many of the foreground pixels are correctly identified, Specificity equals true negative rate and indicates how many of the background pixels are correctly identified, and Precision indicates how many of the predicted foreground pixels are actually foreground pixels. 3.4. Results on phantom Micro-CT data 3.4.1. Qualitative results using 3D reconstructions The qualitative assessment of the segmentation performance using 3D reconstructions, highlights critical elements, such as the importance of small objects and the effect of noisy artifacts on the final result. The 3D reconstructions, illustrated in Fig. 5, were produced using Ray Casting algorithm [46] of VTK 6.3 library in Python 2.7 [47]. It is clear that the proposed algorithm provides superior performance compared to the rest of the methods and finds all the visible objects, providing very similar result to that of the OP. AT is slightly worse than the proposed scheme, missing completely the smallest wire. Otsu and RG provide low quality results with a lot of noise, and AC algorithms are practically nonapplicable, segmenting usually most of the dark background instead of the objects inside it.

3.3. Quantitative evaluation metrics Since the actual foreground pixels constitute less than 4% of the total pixels, the quantitative evaluation of the segmentation results was based on well-established imbalanced evaluation metrics, because typical metrics may not have high enough discrimination ability. These are:

Harmonic F Score ¼ GMean ¼

2⋅Precision⋅Recall Precision þ Recall

3.4.2. Quantitative results The quantitative results of the aforementioned algorithms are presented in Table 2. The proposed IACS method provides superior performance from the rest of the methods under comparison in all metrics, with Harmonic F-Score (HF) (85.1%), G-Mean (86%), Index of Balanced Harmonic F-Score (IBHF) (79.5%) and Index of Balanced G-Mean (IBGMean) (80.3%). Furthermore, IACS is the only method that constantly identifies and segments all the visible objects (14 in total) in the micro-CT OP, even those with marginal size, such as Foil4 and Wire3. Regarding thresholding methods, Otsu thresholding provides mediocre results with HF (47.5%), G-Mean (49.7%), IBHF (53.6%) and IBGMean (76.9%). Its performance is mainly hampered by oversegmentation errors, noise and medium robustness, since it fails completely in some images (Fig. 5(c)). AT has very good performance in all metrics with HF (82.2%), G-Mean (83.4%), IBHF (75.7%) and IBGMean (76.9%). Nevertheless, it does not outperform the proposed method in any of them, providing in general 3–4% lower scores. Its errors are mainly due to noisy artifacts, which appear in a degree analogous to the actual foreground size. Thus, in real images where foreground is larger, this behavior is very likely to lead to oversegmentation and noise, similarly to Otsu. Another very problematic aspect of its behavior is that it misses the very dark areas (with intensity very close to zero), which is a critical error since these areas are certainly foreground. Regarding RG, it has very low performance in all metrics with HF (34.9%), G-Mean (39.7%), IBHF (32.6%) and IBG-Mean (37.4%), although it had a very favorable initialization scheme with seeds manually annotated inside each object. The results showed oversegmentation of larger objects, and undersegmentation or complete omission of smaller objects. This behavior persisted regardless of the density of seeds or their positioning inside the objects. Regarding the AC family of algorithms, it is evident that they perform poorly regardless of the selected implementation. Indicative of the situation is the fact that the best results are obtained by GAC with HF (3.4%), G-Mean (11.3%), IBHF (4.0%) and IBG-Mean (13.5%); numbers that are much lower than those of the proposed IACS. Finally, the proposed method was tested for statistical similarity

(4)

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Sensitivity⋅Specificity

(5)

Table 1 Measured objects' thickness in the micro-CT Optimized Phantom dataset for the applied local threshold values and for the originally proposed global fixed threshold 143. Item

Nominal Thickness Values Manually Measured (μm) Threshold Thickness (μm)

Std Thickness (μm)

Std

Foil 1 Foil 2 Foil 3 Foil 4 Wire 1 Wire 2 Wire 3 Wire 4

250 100 50 20 250 125 50 20

Global FT-143

18.8 5.7 14.1 – 20.3 19.6 – –

Sphere 1 Sphere 2 Sphere 3 Sphere 4 Mesh 1 Mesh 2 Mesh 3

256 102 45 41 247 125 48 –

21.7 19.1 14.1 9 26.9 17 5.4 –

1000

137 129 143 manual 164 174 manual Not obvious 135

264 116 45 – 222,5 95,6 – –

935

112 935

112

1000

135

935

112 935

111

1000

135

935

112 935

111

1000

135

935

112 935

111

100 100 100

145 153 168

102 101 102

21 100.9 23.2 94 24.2 79.2

21.4 23 20.6

(6)

363

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

against the results of the rest of the methods using t-test with 95% confidence level and was found to be statistically different in all five cases, with p-values much lower than 106. 3.4.3. Quantitative results using artificial noise The experimental procedure of section 3.4.2 was repeated under the influence of four types of artificial noise to evaluate the robustness of the proposed method: Gaussian (SNR ¼ 0.1), Poisson, Salt & Pepper (S&P) (density ¼ 0.05) and Speckle (variance ¼ 0.04), to check the behavior of the algorithms under different types of distortions. Noise parameter values were chosen in order to provide typical medium noise intensity. The quantitative results for each noise type and algorithm are presented in Table 3. Under Gaussian noise, IACS provides good performance with HF (84.0%), G-Mean (85.1%), IBHF (76.5%) and IBG-Mean (77.3%), followed by AT with HF (53.9%), G-Mean (57.6%), IBHF (44.7%) and IBGMean (47.7%). Nevertheless, the two methods have an obvious difference of about 30% in all metrics. The rest of the methods provide very poor results with scores below 25%. Under Poisson noise, IACS performs similarly to Gaussian noise. This is not true for AT, whose performance deteriorated dramatically scoring at best 26.9% in G-Mean. RG shows a slight increase in HF (22.6%) and G-Mean (33.1%). Otsu and the AC methods still suffer from very low performance. Under S&P noise, the best results are again provided by IACS with HF (58.6%), G-Mean (59.2%), IBHF (36.7%) and IBG-Mean (37.1%). AT and Otsu provide the second best results, but their scores still remain very low. Indicatively, the highest score for AT is G-Mean (22.6%) and for Otsu is G-Mean (21.5%). The RG and AC methods, provide again extremely low results. Under Speckle noise, the previous pattern is repeated, with IACS scoring HF (76.5%), G-Mean (77.9%), IBHF (54.5%) and IBG-Mean (55.6%), followed in the second place by RG with HF (22.8%), G-Mean (33.5%), IBHF (3.4%) and IBG-Mean (5.1%), in a more than evident difference in performance between the two. The rest of the methods provide even lower results in all metrics, usually below or close to 10%. 3.4.4. Quantitative results of automatic ROI extraction with other segmentation methods The segmentation process was repeated using the intermediate result of the Automatic ROI extraction step but combined with a different segmentation method, namely Otsu, AT and the original CV, in order to investigate both the Automatic ROI Extraction step's influence and the AutoCoefficient CV modification's influence on the final segmentation performance. The results are illustrated on Table 4: The use of the Automatic ROI Extraction step significantly enhances both the Otsu and the original CV algorithms. More specifically, for Otsu, HF has increased by 23.1%, G-Mean by 24.5%, IBHF by 16.9% and IBGMean by 18.0%, while for CV the improvements are much more spectacular with an increase by 74.3% in HF, 70.3% in G-Mean, 73.8% in IBHF and 69.5% in IBG-Mean. However, their performance still remains clearly inferior to that of the proposed IACS in all metrics, a fact that also illustrates the positive impact of the proposed AutoCoefficient modification on the CV algorithm. Finally, the Automatic ROI Extraction step has no impact on the performance of the AT algorithm, which most likely already uses intrinsically the local information through its kernel functionality. 3.5. Results on real Micro-CT data 3.5.1. Qualitative results on 2D image case studies Three 2D image cases were selected in order to illustrate the behavior of IACS, AT and Otsu on real data, shown in Fig. 6 (a)-(c). The cases were chosen in order to illustrate algorithms' behavior on images with different foreground attributes, as it will be explained specifically for each case. For all the algorithms, the same configurations (parameters' values, initialization scheme etc.) with the micro-CT Phantom dataset are used, since that can be considered optimal for micro-CT data.

Fig. 5. 3D reconstructions of the results of each algorithm on micro-CT Phantom dataset. (a) Optimized Phantom, (b) IACS, (c) Otsu thresholding, (d) Adaptive Thresholding, (e) Region Growing, (f) Chan-Vese AC, (g) Geodesic AC and (h) ALR-CV AC.

364

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

Table 2 Quantitative results on micro-CT Phantom. Method

HF

G-Mean

IBHF

IBG-Mean

Objects

IACS Otsu AT RG CV GAC ALR-CV

0.851 0.475 0.822 0.349 0.026 0.034 0.022

0.860 0.497 0.834 0.397 0.090 0.113 0.073

0.795 0.536 0.757 0.326 0.028 0.040 0.020

0.803 0.561 0.769 0.374 0.095 0.135 0.066

14 12 13 12 – 1 –

The best results of each metric are presented in bold. Table 3 Quantitative results on micro-CT Phantom using various types of artificial noise. Noise Type

Method

HF

G-Mean

IBHF

IBG-Mean

Gaussian

IACS

0.840

0.851

0.765

0.773

Otsu AT RG CV GAC ALR-CV

0.083 0.539 0.165 0.027 0.032 0.022

0.148 0.576 0.243 0.092 0.111 0.072

0.105 0.447 0.019 0.029 0.038 0.020

0.186 0.477 0.029 0.097 0.133 0.065

IACS

0.844

0.853

0.767

0.776

Otsu AT RG CV GAC ALR-CV

0.084 0.177 0.226 0.028 0.033 0.022

0.148 0.269 0.331 0.093 0.112 0.072

0.106 0.154 0.035 0.030 0.039 0.020

0.186 0.233 0.051 0.099 0.133 0.065

IACS

0.586

0.592

0.367

0.371

Otsu AT RG CV GAC ALR-CV

0.138 0.136 0.066 0.031 0.025 0.031

0.215 0.226 0.075 0.108 0.087 0.097

0.095 0.111 0.009 0.038 0.026 0.031

0.148 0.185 0.010 0.131 0.087 0.094

IACS

0.765

0.779

0.545

0.556

Otsu AT RG CV GAC ALR-CV

0.019 0.019 0.228 0.034 0.032 0.023

0.085 0.073 0.335 0.113 0.111 0.074

0.025 0.018 0.034 0.041 0.038 0.021

0.112 0.072 0.051 0.136 0.132 0.066

Poisson

S&P

Speckle

Metric

The best results of each metric are presented in bold. Table 4 Quantitative results on micro-CT Phantom using combinations of the Automatic ROI Extraction step with various segmentation algorithms. Method

HF

G-Mean

IBHF

IBG-Mean

IACS Otsu AT Original CV

0.851 0.706 0.822 0.769

0.860 0.742 0.834 0.793

0.795 0.705 0.757 0.766

0.803 0.741 0.769 0.790

The best results of each metric are presented in bold.

conditions. IACS (Fig. 6 (e)) provides accurate and robust performance to all types of foreground. AT (Fig. 6 (h)) introduces a lot of noise, especially in the denser regions, and completely misses some of the darkest areas in the upper left and lower right parts of the image. This is a very critical mistake, since these regions are in reality dense foreground. Otsu (Fig. 6 (k)) introduces some noise, mainly visible in the periphery of the bone, while returning coarser foreground which is visible at larger regions. Case 3 (Fig. 6 (c)) is comprised of similar structures of medium thickness, showing a typical healthy dense bone. IACS (Fig. 6 (f)) performs as previously, returning accurate and robust results, without introducing noise or coarser contours. AT (Fig. 6 (i)) suffers from heavy noise and Otsu method (Fig. 6 (l)) is stable but returns coarser regions

Case 1 (Fig. 6 (a)) is characterized by small and multiple foreground regions, containing many small trabeculae. The image is typical of weak or broken bone. The proposed IACS (Fig. 6 (d)) manages to find most of them, performing robustly without introducing any noise. AT in (Fig. 6 (g)) returns more foreground but with the addition of noisy artifacts results especially in the upper and left part of the image. Otsu (Fig. 6 (j)) seems to provide results somewhere between the two aforementioned methods, being more stable and less noisy than AT, but returning coarser foregrounds than those of the proposed method. Case 2 (Fig. 6 (b)) is showing a combination of dense trabecular bone regions and thin regions and small trabeculae (lower left part), aiming at illustrating the ability to adapt and perform well under mixed foreground 365

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

Fig. 6. Three cases illustrating the behavior of the top three algorithms. (a)–(c) The original images of Case 1, Case 2 and Case 3 respectively. (d)–(f) Results of IACS on the three cases. (g)–(i) Results of Adaptive Thresholding on the three cases. (j)–(l) Results of Otsu on the three cases.

(f). The same overall pattern demonstrated previously on the Phantom and the real 2D Cases is found here as well. Specifically, AT returns a lot of small noisy artifacts. Some of them are clearly visible around the sample and others are inside the spongy bone, not easily visible from a macroscopic view. Otsu, is the worst of three methods, characterized by a big miss-segmentation around the broken area and coarser contours (most clearly visible on the bottom of the sample). On the other hand, IACS provides both clear and accurate results, with tighter contours.

and some small artifacts in the periphery of the bone. In all of the above cases, a similar pattern has been observed where the proposed IACS provides high object identification ability, returning tight and accurate regions without introducing artifacts and noise. AT returns slightly more foreground compared to the rest of the algorithms, but introduces heavy noise. Otsu, introduces some noise and returns most of the foreground but with coarser regions, a deviation which may have severe impact on the calculation of the bone's morphometrical parameters [48].

4. Discussion 3.5.2. Qualitative results using 3D reconstructions In Fig. 7, 3D reconstructions of a real bone sample before (left column) and after compression (right column) are illustrated, using the topthree performing methods, namely IACS (a)-(b), AT (c)-(d) and Otsu (e)-

4.1. Summary of the algorithmic evaluation The proposed IACS framework, provides some important advantages

366

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

over the currently established methods of micro-CT bone segmentation. Specifically, it has enhanced object identification ability even for very small structures as it has been demonstrated on both the micro-CT Phantom dataset and on real data. Secondly, it is an automatic procedure which does not require any real-time user interaction with the data. Although it has three tunable parameters, these are applicationdependent instead of dataset-dependent, and should be changed only once in the case of different image types/applications to adapt the algorithm to the new image modality. As a result, for micro-CT images no parameter selection or optimization is required. Third, the proposed method is highly robust, providing the intended performance regardless of dataset's or slice's details, without introducing spurious artifacts and noise, as it has been illustrated on both Phantom and real cases. This attribute is also maintained to a high degree under all types of noise, where it constantly provides the best results. This consistency is not apparent to the rest of the methods, such as AT that is its main competitor. Fourth, all the above attributes are complemented with high area accuracy. Finally, IACS constitutes an effective implementation of AC algorithms on bone micro-CT data, since typical AC algorithms provide very poor results both qualitatively and quantitatively, as it has been shown in Sec. 3. All the aforementioned attributes represent a significant improvement in the overall current state of micro-CT bone segmentation. The IACS framework has three parameters, two in the Automatic ROI Extraction module, namely BE step's iterations and image patch size, and one in the IAC Evolution module, namely the regularization coefficient of the AC (typical AC parameter). The number of BE iterations depends on the number of backgrounds present in the image. Larger numbers increase the intensity of this processing but may lead to wrong results and small numbers may not remove enough background. Consequently, the typical values should be small positive integers. The patch size is mainly dependent on the image dimensions. An empirical rule for determining the range of appropriate values is logðPatch sizeÞ=logðMax Image dimensionÞ ¼ 0:5 ± 0:1, where Max Image dimension is the largest dimension of the image. For example, in the given micro-CT image dataset the largest image dimension is 1024, so the range of suitable values is between 16 and 64. In Fig. 8, a graph of the F-score and G-mean for the various patch sizes on the micro-CT Phantom dataset is presented, where the optimal value for the patch size is found to be 32. This value satisfies the aforementioned rule. Experimentation for the optimal value maybe required for other image modalities. In the IAC Evolution, the number of AC parameters is reduced from three to one due to the replacement of the two force coefficients with the proposed fully automatic modification of the AC data term. For the remaining regularization coefficient parameter, its values are decimals smaller than one, which are typical values for the Chan-Vese AC [35]. Finally, the independence attribute of the AC algorithms is implemented without the need for any new special parameters. AT is a thresholding method that provides competitive performance on the Phantom dataset, but not as good on the real data. Its main advantages are the semi-automatic nature of the algorithm and the use of local information. AT shows sensitivity to both its quantitative parameters, but the local threshold parameter appears to be the most critical. Regarding the mean intensity calculation method, the algorithm performed well only using actual mean. The optimal values were found by experimentation on the Phantom dataset that was based on changing the value of one parameter while keeping the others constant until the optimal results were found, and then repeating for the rest of the parameters. The results showed that local threshold parameter values close to 21 provide maximum object identification ability but also lead to a lot of noise, and values close to 43 do the opposite. Consequently, the chosen value is 31, striking a good balance between the two extremes. Although AT has high object identification ability, it also introduces a lot of noise, so this benefit cannot be straightforwardly harnessed. Another major drawback is that it misses very dark areas, which are certainly foreground, when any of the above configurations is used. So there is a clear tradeoff between

Fig. 7. 3D reconstruction of a real bone dataset. The left column contains the sample before compression and the right column the sample after compression. (a)–(b) IACS before and after compression. (c)–(d) AT before and after compression (e)–(f) Otsu before and after compression.

367

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

Fig. 8. Values of F-Score and G-mean metrics for different patch sizes on the Phantom dataset.

IACS does not provide as fast execution times as Otsu or AT do, needing about 10 min for the micro-CT Phantom and about 160 min for a real dataset. However, these times are typical for all the applied AC algorithms. Moreover, the experiments were executed using sequential processing and a Matlab implementation. The increase in execution efficiency through the use of a performance oriented language, such as C/ Cþþ, and the use of parallel processing schemes, such as GPU processing, that apply IACS simultaneously on multiple images, will decrease dramatically IACS0 execution times. Such implementations however, deviate from the scope of the current work. Finally, IACS may need parameter tuning in order to be applicable to different image modalities. The Phantom dataset is composed by objects which are a simplification of the actual trabecular bone structures, thus present different characteristics. However, other available phantoms, like densitometric micro-CT phantoms, do not guarantee the dimensions of the geometry and therefore do not allow the evaluation of a clear segmentation error. Applying a simulation of trabecular bone would limit the study to the use of a synthetic phantom with added noise, since a real acquisition wouldn't be available. In the original publication of micro-CT Phantom dataset [2], the authors proposed a fixed threshold of 143 (FT-143) for optimal segmentation, calculated through a procedure that uses the apriori known morphometrical parameters' values. The aluminum objects included in the Phantom have very similar attenuation coefficient to that of the bone and therefore, at the moment, this phantom is still the best option for the selection of the best FT. Any subjective expert based selection will lead to bigger error. Nevertheless, this segmentation technique has intrinsic limitations. For example, the application of the FT-143 on the Phantom completely misses the three smallest visible objects, while IACS misses only one, which is an important advancement on its own. Regarding future directions, IACS0 direct 3D implementation is the first priority that may lead in a decrease in execution times and increase in the accuracy through the use of intra-slide information. Another aim is the application of the current method on more medical image modalities, such as MRI. Thirdly, a more efficient implementation using performance-oriented programming languages and GPU processing techniques is needed, in order to increase execution speed. Finally, since the two modules of the proposed framework are independent, execution of the IAC evolution with manual ROIs or use of other segmentation algorithms with the Automatic ROI Extraction module, should provide interesting future results and applications.

small trabeculae identification and noise introduction. Moreover, its parameters must be manually tuned for each dataset and a form of denoising should be applied, with the risk of losing some of the correctly segmented small trabeculae in the end. Concluding, although AT appears to be a very good candidate for bone micro-CT automatic segmentation, the abovementioned facts reduce its general applicability. Otsu is one of the most popular bone micro-CT automatic segmentation methods. However, from the results in Sec. 3 it appears that it provides medium accuracy and low robustness, due to oversegmentation, introduction of noise and instability even in neighboring slices. Also, its performance deteriorates dramatically under the influence of synthetic noise. As a result, it is not very reliable for large scale applications. In the same manner, RG suffers from undersegmentation and unreliable performance even in neighboring slices, as it was shown in both quantitative and qualitative results. Experimentation with its parameter's values does not provide significant benefits. Moreover, it needs intelligent seed initialization, which is usually done manually and is neither trivial nor fast. In conclusion both of these methods provide medium to low performance under normal or noisy conditions. Typical AC methods provided inadequate results both qualitatively and quantitatively on the Phantom dataset. As a result, these methods cannot be used practically in their original form. However, both quantitative and qualitative results prove that the proposed IACS represents an effective novel algorithm (based on AC algorithmic techniques) for trabecular bone micro-CT data segmentation. 4.2. Application of IACS to other image modalities In order to investigate the general applicability of the proposed IACS, the algorithm was applied for the automatic segmentation of various anatomical regions from different imaging modalities (Fig. 9). In Fig. 9 (a), the automatic ROI extraction result for a CT Axial spine with dimensions 512  512 is presented, where the bone (white) is segmented from the rest of the tissues with very good accuracy (Fig. 9 (b)). The complement of the image was processed in order to have black foreground instead of white (similarly to micro-CT data). In Fig. 9 (c), the ROI of a histological blood cell image with dimensions 500  300 is presented and the corresponding segmentation results in Fig. 9 (d). For both cases, the BE iterations are reduced to 2 and the patch size to 16. The results show high ROI extraction ability, since all the cells are identified and segmented correctly. It should be noted that the aforementioned empirical rule for determining the suitable patch size value in Section 4.1 is also verified in both images from different modalities. The above cases, serve as a brief demonstration of the capabilities of the IACS framework on different medical image modalities. By tuning two intuitive parameters within a very limited range, the results can be considered more than satisfactory and provide a sound argument in favor of the general applicability of the proposed IACS scheme.

5. Conclusion In this paper a novel segmentation framework for bone micro-CT images was proposed. The framework consisted of two independent modules: Automatic ROI Extraction and Independent AC Evolution. The proposed IACS framework was tested against several established segmentation methods namely, Adaptive Thresholding, Otsu Thresholding, Region Growing, Chan-Vese AC, Geodesic AC and ALR-CV AC. The method was evaluated quantitatively on a micro-CT Optimized Phantom

4.3. Limitations and future directions The study also presents some limitations. Specifically, the proposed 368

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370

Fig. 9. Application of IACS in different medical image modalities. (a) ROI of CT-spine image and (b) the segmentation result (c) ROI of the histological blood cell image and (d) the segmentation result. € [12] C. Ohman, et al., Mechanical testing of cancellous bone from the femoral head: experimental errors due to off-axis measurements, J. Biomech 40 (11) (Jan. 2007) 2426–2433. [13] M.L. Bouxsein, et al., Guidelines for assessment of bone microstructure in rodents using micro–computed tomography, J. Bone Min. Res 25 (7) (Jul. 2010) 1468–1486. [14] S.-W. Lee, et al., Stem cell-mediated accelerated bone healing observed with in vivo molecular and small animal imaging technologies in a model of skeletal injury, J. Orthop. Res. 27 (3) (Mar. 2009) 295–302. [15] M. Ding, et al., Accuracy of cancellous bone volume fraction measured by micro-CT scanning, J. Biomechanics 32 (3) (1999) 323–326. [16] I. Lima, et al., Ethanol bone evaluation using 3D microtomography, Micron 39 (5) (Jul. 2008) 617–622. [17] J.-T. Hsu, et al., The assessment of trabecular bone parameters and cortical bone strength: a comparison of micro-CT and dental cone-beam CT, J. Biomech. 46 (15) (Oct. 2013) 2611–2618. [18] G. Michel, et al., Micro-ct analysis of radiation-induced osteopenia and bone hypovascularization in rat, Calc. Tiss. Int. 97 (1) (Jul. 2015) 62–68. [19] S. Tassani, G.K. Matsopoulos, The micro-structure of bone trabecular fracture: an inter-site study, Bone 60 (Mar. 2014) 78–86. [20] L. De Schaepdrijver, et al., In vivo longitudinal micro-CT study of bent long limb bones in rat offspring, Repr. Tox. 46 (Jul. 2014) 91–97. [21] T. Hara, et al., The influence of microcomputed tomography threshold variations on the assessment of structural and mechanical trabecular bone properties, Bone 31 (1) (Jul. 2002) 107–109. [22] A. Otsu, Threshold selection method from gray-level histograms, IEEE Trans. Syst. Man Cybern. 9 (1) (Jan. 1979) 62–66. [23] P. Dong, et al., 3D osteocyte lacunar morphometric properties and distributions in human femoral cortical bone using synchrotron radiation micro-CT images, Bone 60 (1) (Mar. 2014) 172–185. [24] H.R. Buie, et al., Automatic segmentation of cortical and trabecular compartments based on a dual threshold technique for in vivo micro-CT bone analysis, Bone 41 (4) (Oct. 2007) 505–515. [25] M. Sezgin, B. l Sankur, Survey over image thresholding techniques and quantitative performance evaluation, J. Electron. Imaging 13 (1) (2004) 146–168. [26] W.K. Pratt, Digital image processing, fourth ed., John Wiley & Sons, Los Altos, California: USA, 2007. [27] H. Scherf, R. Tilgner, A new high-resolution computed tomography (CT) segmentation method for trabecular bone architectural analysis, Am. J. Phys. Anthropol. 140 (1) (2009) 39–51. [28] A.C. Jones, et al., Analysis of 3D bone ingrowth into polymer scaffolds via microcomputed tomography imaging, Biomat. 25 (20) (Sept. 2004) 4947–4954. [29] S. Tassani, et al., Influence of segmentation on micro-CT images of trabecular bone, J. Microsc. 256 (2) (Nov. 2014) 75–81. [30] M. Kass, et al., Snakes: active contour models, Int. J. Comp. Vis. 1 (4) (Nov. 1988) 321–331. [31] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comp. Phys. 79 (1) (Nov. 1988) 12–49. [32] D. Li, et al., Active contours driven by local and global probability distributions, J. Vis. Comm. Im. Repres. 24 (5) (Jul. 2013) 522–533.

dataset, under normal conditions and under four types of artificially introduced noise and was shown to provide superior performance in terms of accuracy, robustness and object identification ability. Moreover, qualitative evaluation was made using 3D reconstructions of the microCT Phantom and of a real case, and using three 2D real data case studies. Finally, in order to demonstrate its general applicability, IACS was applied in two other medical image modalities (CT, Histology) providing more than acceptable results without any significant changes in its configuration. Acknowledgments We would like to thank Prof. Fabio Baruffaldi, from the Laboratorio di Tecnologia Medica, Istituto Ortopedico Rizzoli, Bologna, Italy for providing the micro-CT Phantom and the real trabecular bone images. References [1] D.P. Clark, C.T. Badea, Micro-CT of rodents: state-of-the-art and future perspectives, Phys. Med. 30 (6) (Sep. 2014) 619–634. [2] E. Perilli, et al., A physical phantom for the calibration of three-dimensional X-ray microtomography examination, J. Microsc. 222 (2) (May 2006) 124–134. [3] E. Perilli, et al., MicroCT examination of human bone specimens: effects of polymethylmethacrylate embedding on structural parameters, J. Microsc. 225 (2) (Feb. 2007) 192–200. [4] D.D. McErlain, et al., Study of subchondral bone adaptations in a rodent surgical model of OA using in vivo micro-computed tomography, Osteoarth. Cart. 16 (4) (Sep. 2008) 458–469. [5] A. Laib, et al., The temporal changes of trabecular architecture in ovariectomized rats assessed by MicroCT, Osteop. Int. 12 (11) (Nov. 2001) 936–941. [6] C. Cowan, et al., MicroCT evaluation of three-dimensional mineralization in response to BMP-2 doses in vitro and in critical sized rat calvarial defects, Tiss. Engin. 13 (3) (Mar. 2007). [7] S.T. Ho, D.W. Hutmacher, A comparison of micro CT with other techniques used in the characterization of scaffolds, Biomat 27 (8) (Mar. 2006) 1362–1376. [8] D.W. Dempster, et al., Effects of daily treatment with parathyroid hormone on bone microarchitecture and turnover in patients with osteoporosis: a paired biopsy study*, J. Bone Min. Res. 16 (10) (Oct. 2001) 1846–1853. [9] J. Waarsing, et al., Longitudinal micro-CT scans to evaluate bone architecture, J. Musculoskelet. Neuronal Interact. 5 (4) (Oct.-Nov. 2005) 310–312. [10] J.H. Kinney, et al., In vivo, three-dimensional microscopy of trabecular bone, J. Bone Mins Res. 10 (2) (Feb. 1995) 264–270. [11] P. Rüegsegger, et al., A microtomographic system for the nondestructive evaluation of bone architecture, Calc. Tiss. Int. 58 (1) (Jan. 1996) 24–29.

369

V.Ch. Korfiatis et al.

Computers in Biology and Medicine 87 (2017) 358–370 [40] P.L. Rosin, Unimodal thresholding, Patt. Rec. 34 (11) (Nov. 2001) 2083–2096. [41] N. Coudray, et al., Robust threshold estimation for images with unimodal histograms, Patt. Rec. Lett. 31 (9) (Jul. 2010) 1010–1019. [42] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Patt. An. Mach. Intell. 12 (7) (Jul. 1990) 629–639. [43] S. Yonggang, W.C. Karl, Real-time tracking using level sets, in, IEEE Comput. Soc. Conf. Comp. Vis. Patt. Rec. 2 (2005) 34–41. [44] Y. Tian, et al., Active contour model combining region and edge information, Mach. Vis. Appl. 24 (1) (Jan. 2013) 47–61. [45] Bruker microCT. Available: http://bruker-microct.com/products/ctan.htm. [46] S.D. Roth, Ray casting for modeling solids, Comput. Graph. Image Process. 18 (2) (Feb. 1982) 109–144. [47] The Visualization Toolkit. Available: http://www.vtk.org/. [48] M. Ding, et al., Age-related variations in the microstructure of human tibial cancellous bone, J. Orthop. Res. 20 (3) (2002) 615–621.

[33] K. Zhang, et al., Active contours with selective local or global segmentation: a new formulation and level set method, Im. Vis. Comp 28 (4) (Apr. 2010) 668–676. [34] V. Caselles, et al., Geodesic active contours, Int. J. Comp. Vis. 22 (1) (Feb. 1997) 61–79. [35] T.F. Chan, L.A. Vese, Active contours without edges, IEEE Trans. Im. Proc. 10 (2) (Aug. 2001) 266–277. [36] A. Tsai, et al., Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification, IEEE Trans. Im. Proc. 10 (8) (Aug. 2001) 1169–1186. [37] V.C. Korfiatis, et al., Automatic local parameterization of the Chan Vese active contour model's force coefficients using edge information, J. Vis. Comm. Im. Repr. 29 (May 2015) 71–78. [38] D. Mumford, J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems, Comm. Pure Appl. Math 42 (5) (Jul. 1989) 577–685. [39] G.W. Zack, et al., Automatic measurement of sister chromatid exchange frequency, J. Histochem. Cytochem 25 (7) (July 1977) 741–753.

370