Quantitative analysis of bone microvasculature in a mouse model using the monogenic signal phase asymmetry and marker-controlled watershed

Thethree-dimensional (3D) imaging and quantitative analysis of bone microvasculature are important to describe angiogenesis involvement in bone metastatic processes. Here, we propose an algorithm based on marker-controlled watershed for the 3D segmentation of vessels and bone in mouse bone imaged with a contrast agent using synchrotron radiation micro-computed tomography (SR-μCT). Markers were generated using hysteresis thresholding and morphological filters, and the control surface was constructed using the monogenic signal phase asymmetry. The accuracy and robustness of the proposed method were evaluated on a series of synthetic volumes generated to mimic the vessel, bone and background structures. Different contrast between different structures, as well as different noise levels were considered. A series of multi-class synthetic volumes were segmented using the proposed method, and the overall segmentation quality was evaluated using the Matthews correlation coefficient (MCC) by comparing to the ground truth. Additionally, we evaluated the segmentation of thin structures under various levels of Gaussian noise. The simulation study indicated that the algorithm was performant in multi-class segmentation with different contrast, noise, and thickness. The algorithm was applied to images of bone from a mouse model of breast cancer bone metastasis acquired using SR-μCT. The segmentation quality was evaluated using the Dice coefficient and the MCC by comparing to manual segmentation. The proposed method performed better than hysteresis thresholding and marker-controlled watershed using the magnitude of the gradient as control surface. Several quantitative parameters on bone and vessels were extracted, including bone volume fraction (BV/TV), vessel volume fraction (VV/TV) and the mean vessel thickness (VTh). The bone volume fraction (BV/TV) was significantly lower in the metastatic group compared to the healthy group. This demonstrated the effectiveness of the algorithm for the study of bone and vessel microstructures in mouse model.


Introduction
Breast cancer is the most frequently diagnosed cancer in women worldwide (Bray et al 2018). In terms of mortality, there are approximately 508 000 deaths annually around the world (Tulotta et al 2019). The major morbidity and mortality are caused by bone metastases, generally leading to osteolytic lesions. Around 70% of breast cancer patients develop to the advanced phase where tumors metastasize to bone (Akhtari et al 2008). Once metastases locate in bone, the existing treatments are not able to cure the disease, and the median survival rate is about 24-65 months (Lote et al 1986, Nutter et al 2014, Catarina et al 2017.
Breast cancer bone metastases can not only bring about bone destructions but also facilitate the formation of undesirable vascularization (Nyangoga et al 2011). To investigate and have a better understanding of the pathology, the three-dimensional (3D) imaging and analysis of bone and vasculatures are necessary. 3D x-ray micro-computed tomography (μCT) is a tool of choice to image trabecular bone microarchitecture (Kuhn et al 1990) as well as cortical bone and its Haversian network (Cooper et al 2003, Bousson et al 2004. In the study of assessing the effects of estrogen deficiency on bone, 3D x-ray μCT was used to image and analyze the rat bone structures (Sharma et al 2018). Nevertheless, the protocol did not image the blood vessel to assess its alterations caused by the estrogen deficiency. In addition, 3D x-ray μCT can be used to visualize the vascular architecture with contrast agent (Moore et al 2003, Zhang et al 2005. However, these studies require to decalcify the bone. Therefore, μCT is only used to image and analyze either bone microstructure or vascular architecture separately, and these studies do not permit a complete analysis of the relationships between the bone and vasculature. Synchrotron radiation μCT (SR-μCT) has significant advantage of acquiring high spatial resolution images with a high signal-to-noise ratio, compared to the standard μCT, due to the high photon flux of the synchrotron source (Salomé et al 1999). SR-μCT has been applied to the quantitative analysis of human trabecular bone up to the micrometer scale (Peyrin et al 1998, Larrue et al 2011. In addition, SR-μCT coupled with the use of a contrast agent permitted to visualize simultaneously the 3D bone microstructures and vascular networks in mouse (Schneider et al 2009) or rat (Prisby et al 2011). Simultaneous visualizations of the calcified bone microstructure and blood vessels were also demonstrated using phase contrast synchrotron radiation computed tomography (SR-PCT) without any contrast agent in mice model (Núñez et al 2017). However, in this study, the field of view (FOV) was limited due to the small pixel size used, and vessels were segmented by hand due to the relatively weak phase contrast in the vessels limiting the applications of the protocol to large datasets.
The simultaneous imaging of bone and vessel is often followed by performing segmentation to separate bone and vessel compartments. The bone and vessel structures of mice, simultaneously imaged using SR-μCT with a contrast agent, have been segmented automatically using the global thresholding (Schneider et al 2009). We have previously proposed 3D region growing to segment rat bone vascularization (Langer et al 2010, Prisby et al 2011 in SR-μCT images associated to a contrast agent. Specifically, vessels were segmented using hysteresis thresholding (a particular case of region growing with two thresholds) to overcome partial volume effects.
The aim of this work is to transfer the previous protocol to mice in order to make more models of pathologies available. However, this transition is not straightforward. Firstly, as opposed to in rat bone, vessels may appear to be in contact with the bone surface in mice, precluding the correct segmentation of bone and vessels using the previously proposed hysteresis thresholding based protocol, at 3.5 μm resolution.
To address these challenges, we proposed to use the marker-controlled watershed algorithm (Beucher and Meyer 1990) in conjunction with the monogenic signal phase asymmetry to segment and quantify the bone microstructure and vascular network in SR-μCT images of mice. The proposed method was not only assessed on synthetic volumes, but also validated on experimental dataset by comparing to the manual segmentation. Segmentation using the proposed method was quantitatively compared to the hysteresis thresholding based method and gradient based marker-controlled watershed. To quantitatively analyze vessels and bone in the healthy and metastatic groups, several parameters were extracted to characterize bone microstructures and vasculatures.

Image segmentation
In this study, the aim of the segmentation method is to classify voxels into three classes: bone, vessels and background.
The watershed algorithm can be used to separate vessels appearing to be touching the bone surface, without the need for post-segmentation merging of regions. The classical watershed transform is often used to segment an image into multiple objects by detecting the 'catchment basins' and 'ridge lines'. In this method, the 'flooding' begins from the regional minimums on the 'control surface', which is treated as a topographic map and can be the magnitude of gradient. The ridge lines present in the image will then separate catchment basins when the different sources of water are merging. However, in most cases, there are many undesired regional minimums due to noise or natural variations, which lead to over-segmentation (Bhabatosh 2011), requiring a posteriori merging of the segments. One way to avoid this problem is the marker-controlled watershed method, which initializes the watershed from the already identified markers (Beucher and Meyer 1990). The method saves running time by reducing the number of iterations and prevents results from over-segmentation by simplifying the merging of the regions.
The proposed protocol was based on the marker-controlled watershed algorithm. In this work, we generated the markers of the different classes by using hysteresis thresholding and morphological filters, and constructed the control surface from the calculation of the monogenic signal phase asymmetry of the original image.

Markers generation
An important step of the marker-controlled watershed algorithm is to generate the marker image. To achieve coverage of all connected components in each class and minimize the false positives, the markers can be created using some feature detection methods or by hand (Wang and Vallotton 2010). In this study, the marker generation procedure is introduced by using an example of mouse bone in SR-μCT image, as shown in figure 1. The markers of the three classes are generated using the hysteresis thresholding and morphological filters. The background and vessel markers are given by the cyan and magenta arrows, respectively. The initial segmentation is obtained using the hysteresis thresholding with two thresholds. The selections of the low and high threshold of hysteresis thresholding depend on the gray level distributions in the SR-μCT image. In this study, hysteresis thresholds were determined manually according to the testing results. We applied the same parameters on all samples, since all images had a similar gray level distribution. Then, the initial segmentation is followed by one iteration of morphological thinning to reduce the false-positive rate. The procedure to obtain bone marker is given by the yellow arrows. The whole foreground structure, containing both bone and vessel, is first segmented using hysteresis thresholding. Then, the dilated vessel markers are subtracted to obtain the initial bone markers. Since the aim of morphological dilation was to remove 'halo' structures caused by subtraction, a´3 3 3 voxels ( m m ḿ10.5 m 10.5 m 10.5 m) of structuring element was enough in this study. Lastly, the isolated small particles are removed using morphological opening, and minimizing false positives is achieved by using morphological thinning. The final marker image is generated by adding the background markers, vessel markers, and bone markers, in cyan, magenta, and yellow, respectively.

Control surface generation
The other essential step of marker-controlled watershed algorithm is the generation of a control surface from the original image. Classically, the magnitude of the gradient is used as the control surface, which is related to intensity changes in the original image (Bhabatosh 2011). However, in our case, the contrast at the interfaces between bone and vessels is relatively weak compared to the bone-background and the vessel-background interfaces, leading to weaker gradient magnitude at these edges. To alleviate this problem, we proposed the local phase asymmetry of the 3D monogenic signal as the control surface to improve the edge detection at the low contrast interfaces. The phase asymmetry measurement can be used to quantify whether the signal is locally edge like or line like. Since we aim to detect the interfaces between different structures, we use this property to select only the edge like structures for the control surface.
The 3D monogenic signal can be seen as the isotropic and multidimensional extension of the 1D analytic signal (Felsberg and Sommer 2001), and can be defined by the combination of the original signal and its three Riesz transform components (Chenouard and Unser 2012): 3 is the spatial coordinate. R stands for the Riesz transform, and its operator component is characterized by its frequency response: [ where,  represents the Fourier transform operator, j is the imaginary unit, 3 is the angular frequency variable conjugate to x. a gives the 3D direction of the Riesz transformed components, corresponding to the directions of the basis vectors. In addition, structure features must be selected with an appropriate scale using a band-pass filter, for example a log-Gabor filter with a shape parameter to govern the bandwidth of the passband, since structure in an image is generally scale dependent (Bridge 2017). Thus, the 3D monogenic signal ( ) x f mg is constructed as: stands for log-Gabor filter in spatial domain, and * represents the 3D convolution operation. The 3D monogenic signal can also be constructed as the combination of its even and odd components (Rajpoot et al 2009): where the even component is To detect edges in the original image, the multiscale phase asymmetry is measured by Bridge (2017): where ⌊.⌋ denotes an operator replacing negative values with zero, and { } l i are a set of center-wavelengths ( ) / l p w = 2 i i of the log-Gabor filters to define the scales. T is a threshold to suppress noise,  is small number to avoid division by zero. In this study, we set multiscale values of l = 6, 7, 8, i = T 0.18 and =  0.001.

Segmentation evaluation
To evaluate segmentation quality, the Dice coefficient and Matthews correlation coefficient (MCC) were measured by comparing to manual segmentation or ground truth.

Dice coefficient
The Dice coefficient (also known as Dice similarity coefficient and F1 score) is a statistic parameter to assess spatial overlapping between the segmentation and ground truth (Zou et al 2004). Its expression is given by: where | | A and | | B are the cardinalities of the two sets A and B, representing respectively the segmentation and ground truth (Aaron et al 2020).

Matthews correlation coefficient
The MCC (Matthews 1975) is widely used as a measure of quality of classifications to indicate how well observations (ground truth) and predictions (segmentation) agree.
= MCC 1 indicates a perfect classification, = MCC 0 indicates no better than random prediction, and = -MCC 1 indicates a total misclassification (Matthews 1975). In its original form, MCC applies to measure the quality of two-class classifications. However, a multi-class MCC has been proposed (Gorodkin 2004): where ijth entry C ij is the number of elements of observation class i assigned to class j (Jurman and Furlanello 2010): For instance, in a 3D image, C kk indicates how many voxels are correctly segmented in class k.

Numerical simulations
In order to examine the accuracy and robustness of the proposed method, as well as its generalization to other applications, four kinds of 3D synthetic models were created to simulate real datasets (figure 2). Each volume was an 8 bits 3D image with 256×256×150 voxels. The size of synthetic phantom was determined by striking a balance between the size of real bone in SR-μCT images and computational cost. The gray levels of compartments in synthetic phantoms were set based on the real SR-μCT images. For instance, vessel and background were set at 255 and 0 in phantoms respectively, due to the contrasted vessel (bright) and background (dark) in the real SR-μCT images. As for the gray level of bone, it was set at different values between 0 and 255 to mimic the various contrast. The virtual phantoms were created based on geometrical cylinders and hourglass solids assigned 0, 128, and 255 greyscale values in their initial instances ( figure 2). Additionally, in model-1, model-2, and model-3, the diameters of the widest and narrowest hourglass solids were assigned 34 and 1 voxels, respectively. In model-4, the narrowest gaps between the bright and gray structures were set to 1 voxel. Specifically, model-1 was designed to simulate 3D vessel structure (gray level 255) and bone porosity (background, gray level 0) with various thicknesses, embedded in the bone (gray level 128), as shown in figure 2(a). Model-2 mimicked that vessels (gray level 255) were in contact with bone surface (gray level 128) and background (gray level 0) at the same time, as shown in figure 2(b). Model-3 was designed to simulate an isolated structure (bone or vessel) in the real-world dataset, as shown in figure 2(c). Model-4 described a situation that a vessel (gray level 255) was passing through a canal (porosity) with varying diameters (gray level 0), as shown in figure 2(d). The bright cylinder in the center simulated vessel, and the surrounding structure was the synthetic bone (gray level 128).
To investigate the impact on the segmentation of different contrast levels between bone and vessels in the acquired images, a series of volumes were generated by varying the relative contrast of bone. In model-1, model-2 and model-4, the gray levels of cylinders, which simulated bone structures in the real data, were set to increase starting from the intensity of 48-208 by step of 10. In model-3, the gray level of the hourglass structure was increasing from the intensity of 48-208 by step of 10, and 255.
For the purpose of mimicking the partial volume effect, which is a reduction of gray level in small object and interface between two classes, the phantom images were downscaled and re-upscaled by a factor 4 with the bilinear resampling method. The synthetic phantoms without partial volume effects were used as reference images (ground truth).
In addition, to explore the impact of noise, different levels of Gaussian noise were added to each volume. The noise level was defined by the standard deviation of Gaussian noise. Noise was added with levels in a range of 10-70 with increments of 10, producing 63 instances of model-1, model-2, and model-4, and 70 instances of model-3.
2.4. Application to experimental data 2.4.1. Sample preparation To validate the proposed algorithm on experimental data, we selected two groups of healthy and metastatic mice bone (10 samples per group) from a project investigating the effects of anti-angiogenic drugs (vatalanib: 100 mg kg −1 orally daily; bevacizumab: 5 mg kg −1 intraperitoneally twice a week) on metastasis formation. Specifically, 8 week-old female Balb/c nude mice were injected with luciferase-expressing human B02 breast cancer tumor cells (Fradet et al 2011), then the preventive treatments were performed. In this model, animals usually develop bone metastases 18 d after the injection of tumor cell (Bachelier et al 2014). Thus, the mice were sacrificed on 22 d after the injection. Next, all mice were injected with a contrast agent (barium sulfate) for the vascular imaging. Afterwards, mice tibiae were dissected, embedded in PMMA, and kept for SR-μCT imaging.

Image acquisition
To image the bone microstructure and vasculature simultaneously in 3D, we used contrast agent SR-μCT. Imaging experiments were performed at the European Synchrotron Radiation Facility, Grenoble, France, on the ID19 beamline. The schematic of the SR-μCT imaging setup is shown as figure 3. A 2048×2048 pixel CCDbased detector with effective pixel size of 3.5 μm was used to record images at evenly spaced angles of view over a 360°rotation. The exposure time and x-ray energy were set to 0.15 s and 26 keV, respectively. The acquisition of 2000 radiographs of one sample lasted approximately 8 min. After the image acquisition, 3D 32 bits images of 2000×2000×1200 voxels were reconstructed using a filtered back projection algorithm yielding a cylindrical FOV of 7 mm diameter. The processing involved corrections to take into account beam inhomogeneity, response of the detector and decrease of current during acquisition. From previous experiment, it can be estimated that the actual spatial resolution is of the order of two voxels. Then, the bit depth of reconstructed images were reduced to 8 bits from 32 bits. Finally, the proposed algorithm was applied in the whole SR-μCT images to segment bone and vessels.

Quantitative parameters
To quantitatively analyze the differences between the healthy group and metastatic group, some anatomical parameters were extracted from the data. The total volume (TV) of the sample was defined as the volume inside the outer contour of the bone. TV, bone volume (BV) and vessel volume (VV) were measured by counting voxels in the corresponding compartments. In order to compare the volume-based fractions in the healthy and metastatic groups, the normalized ratios BV/TV and VV/TV were calculated. In addition, the local vessel thickness was measured as the diameter of the largest sphere fitting in the vessel at any point of its length (Hildebrand andRüegsegger 1997, Langer et al 2010), and its average VTh was reported.

Statistical analysis
Statistical tests were performed on measured BV/TV, VV/TV and VTh, in the healthy and metastatic groups, separately. Lilliefors test was used to verify the normal distributions of parameters in each group. In a case of normal distribution, ANOVA F-test was used to test if there was a significant difference between the two groups. We considered the difference statistically significant at p-value<0.05 level.

Results
In this section, we first examine the proposed method of phase asymmetry based marker-controlled watershed on the 3D synthetic volumes. Then we apply the proposed method to the real SR-μCT dataset of mouse bone, and extract the quantitative parameters to characterize bone and vessel structures.

Assessment of the method on synthetic data
We first report two examples of the segmentation of model-1 obtained with various contrast and additive Gaussian noise levels, as shown in figures 4, 2D slices from the 3D volumes. Qualitatively, the generated markers and control surface under a relatively low contrast and high noise level were less robust as figure 4(g), leading to a worse watershed segmentation as figure 4(h).
Specifically, to explore the impact of image contrast and the level of noise on segmentation, we generated 189 synthetic volumes for three-class models (model-1, model-2 and model-4) and 70 volumes for a two-class model (model-3), as explained above. Multi-class structures of each synthetic volume were segmented using the proposed method. In figure 5 we show nine examples of segmentation for each model, obtained under various contrast and noise levels. The panels delineate a representative volume of interest (VOI) for each phantom, obtained at contrast levels of 48, 128 and 208 (or 255) and at noise levels of 10, 30 and 70. Where the algorithm failed to yield a segmentation where the structure was recognizable, which happened only at rather extreme noise levels, we report 'no result'. Qualitatively, the quality of segmentation decreased with the increase of noise level and the decrease of contrast.
Quantitatively, the MCC value for each multi-class segmentation was calculated according to equation (10) as shown in figure 6, which illustrates the relationships between MCC and noise levels, separated according to phantom type and contrast. MCC values were plotted over a range of noise levels at each contrast, excluding extremely low values obtained from the failed segmentation as before ('no result', figure 5). For example, in figure 6(a), when the noise level is higher than 30, with a bone contrast between 48 and 208, the segmentation algorithm seems to reach its limit. This noise level is extreme compared to the noise in the experimental data, however. In summary, for three-class models of figures 6(a), (b) and (d), in a fixed noise level, the highest MCC values were obtained at contrasts of 128 and 148. This indicates that the proposed protocol works best when the intermediary class has gray levels close to half of the maximum gray level. This can serve as a guideline for optimizing the imaging conditions. For the two-class model of figure 6(c), MCC reaches the highest value at the highest contrast of 255, as expected. Figure 7 illustrates the capability of thin structure segmentation under a various levels of Gaussian noise. For the purpose of quantifying the thin structure segmentation, we measured 3D local thickness of structure. Specifically, we selected the hourglass structure from the VOI of ground truth as a foreground to generate a new binary volume, and measured its 3D local thickness map. Secondly, the 3D local thickness map was masked with hourglass structure segmentation of figures 7(a)-(d), and results are shown as figures 7(e)-(h). As expected, the segmentation quality of thin structure was continuously degraded with increasing noise as the disconnection of structure visible. Lastly, the numbers of voxels at different thicknesses, from the figures 7(e)-(h), were plotted in the figure 7(i). In a thickness range of 3-9 voxels (thin structure), it indicates that thin structure can be well segmented at noise levels of less than 20, using the proposed method. For the segmentation of thick structure (thickness of 10 voxels or more), noise levels have little effects on the number of voxels at each thickness. The thick structures can be well segmented even at the noise level of 60 using the proposed method.
3.2. Application to SR-μCT mouse data 3.2.1. Image processing Segmentation results of a representative VOI of an original image are shown in figure 8. A zoom on a 2D slice from the 3D original volume is shown as figure 8(a). The red arrow shows the relatively low contrast interface between bone and vessel. The control surfaces were generated using two 3D edge detection methods: conventional image gradient magnitude in figure 8(b) and phase asymmetry in figure 8(c). As for the conventional gradient, the first-order derivative of an image was approximated by finite differences with convolution kernels. Conventional gradient was not able to properly resolve the low contrast boundary between bone and vessel indicated by the arrow, which was conversely well detected by the phase asymmetry signal. The selection of scales for phase asymmetry depends on image structures. In our case, phase asymmetry with lower scales detects too much edge detail so that useful structures can not be recognized. Conversely, the method with higher scales misses too much detail. As it appears in figure 8(c), even though there are some spurious edges detected in homogeneous areas with the phase asymmetry method, they are always covered by a marker in solid areas. Hence they do not influence the marker-controlled watershed segmentation. Qualitatively, the result of segmentation using the phase asymmetry based marker-controlled watershed was much improved compared to hysteresis thresholding or gradient based watershed segmentation. Hysteresis thresholding requires spatial separation between the different compartments. Since this is not the case here, sufficiently low second thresholds can not be selected to correctly segment the compartments, which leaves a small gap between the structures ( figure 8(g)). For the watershed algorithm, the conventional gradient did not yield sufficient signal at bone to vessel interfaces to properly separate the two compartments ( figure 8(b)).

Segmentation quality evaluation
The segmentation quality was evaluated quantitatively using the Dice coefficient and the MCC by comparing to manual segmentation. A representative volume was manually segmented using VGStudio Max (Volume Graphics GmbH, Heidelberg, Germany). Dice coefficients corresponding to vessel, bone and background compartments were calculated separately, and an overall segmentation quality was also evaluated using multiclass MCC, as shown in the table 1. Especially, there were substantial improvements at the overall segmentation (MCC=0.94) using the proposed phase asymmetry based marker-controlled watershed, comparing to the hysteresis thresholding (MCC=0.77) as well as gradient based marker-controlled watershed (MCC=0.88).

Statistical analysis
To compare the healthy sample with the metastatic sample qualitatively, 3D volume renderings of segmentation using the proposed method are shown as figure 9. The metastatic sample in figure 9(d) shows evident large bone lesions, compared to the healthy sample in figure 9(a). In addition, the vessels in the metastatic sample are more numerous and thicker than those in the healthy sample, as evidenced in the comparison between figures 9(e) and (b). This may indicate an increased and abnormal vascularization due to the bone metastases, as expected. The plots of BV/TV, VV/TV and VTh in the healthy and the metastatic groups are shown as figure 10. Due to the low number of samples, all data points were plotted along with the box plots. Normal distributions of BV/TV, VV/TV and VTh in each group were verified using the Lilliefors test (p-value>0.05). ANOVA F-test (p-value<0.05) was used to test if there is a significant difference between the healthy and the metastatic groups. BV/TV (median: 0.2673) is significantly lower in the metastatic group compared to in the healthy group (median: 0.3263) as figure 10(a) (ANOVA F-test, p-value=0.012), which is consistent with the apparent large metastatic lesions visible in figure 9 (d). In addition, there is no significant difference on parameters VV/TV (p-value=0.715) and VTh (p-value=0.469) as figures 10(b), (c), possibly due to the limited sample size.

Discussion and conclusion
We presented a method to segment bone and vessels using the marker-controlled watershed and the monogenic signal phase asymmetry. The marker-controlled watershed was used to address the problem that vessels may appear to be in contact with the bone surface in mouse model precluding the correct segmentation of bone and vessels using the previously proposed hysteresis thresholding based protocol. In this algorithm, three classes of markers, i.e. bone, vessels and background were generated using the hysteresis thresholding and morphological filters, minimizing the false positives in the connected components of each class. In addition, we proposed the monogenic signal phase asymmetry as the control surface to improve edge detection at the relatively weakly contrasted bone and vessel interfaces.
To examine the accuracy and robustness of the method, as well as its generalization to other applications, we created 3D synthetic models to simulate a real dataset. Since there might be various contrast ratios between vessels and bone, physically came from the different types of contrast agent, x-ray energy, and range for 8 bits resampling in the real experiment, we investigated the impact of image contrast on segmentation quality. In addition, we considered the different noise levels and the segmentation of thin structure. After the multi-class synthetic volumes were segmented using the proposed method, the overall segmentation quality was evaluated using the MCC by comparing to the ground truth. The results indicated that the proposed protocol could work better at the contrast of 128 and 148 in the three-class model, and at the contrast of 255 in the two-class model. In addition, we studied the capability of thin structure segmentation with the proposed method, under various levels of Gaussian noise. The results showed that thin structures (thickness of 3-9 voxels) could be well segmented at noise levels of less than 20, and thick structures (thickness of 10 voxels or more) could be well segmented even at the noise level of 60 using the proposed method. In summary, the proposed segmentation method appears to be robust with respect to noise, contrast and thin structure segmentation. Thus, it is a good candidate for a more general application of multi-class segmentation.
We applied the proposed method to the real dataset of mouse bone with breast cancer bone metastases, acquired using SR-μCT with contrast agent. The segmentation quality was evaluated using the Dice and the MCC by comparing to manual segmentation. The results showed substantial improvements at each single class segmentation (Dice) and overall segmentation (MCC) using the proposed method, compared to the hysteresis thresholding based method as well as to the gradient based marker-controlled watershed. In addition, for the purpose of quantitative analysis in the healthy and metastatic groups, we extracted three parameters (BV/TV, VV/TV and VTh) to characterize bone and vessels structures. We found that the BV/TV was significantly lower in the metastatic group compared to in the healthy group. Thus, the statistical analysis revealed that BV fraction decreased due to the large bone lesions, all in line with previous studies. In addition, there was no significant  difference between the healthy and metastatic groups on parameters of VV/TV and VTh, possibly due to the limited sample size. The application to the real dataset demonstrated the utility of the algorithm for the study of bone and vessel microstructures in mouse models. However, in this study, the statistical analysis might be limited by a small sample size (10 samples per group), leading to a low statistical power and less conclusive results. Further biological study can take this factor into account and add further samples to provide a more statistical basis. The proposed method can be used to gain a better understanding of role of vessels in cancer pathology and can help us better describe angiogenesis involvement in bone metastatic processes. In addition, the proposed protocol can be applied to other possible applications.