Abstract
This paper presents a new strategy for the segmentation of trabecular bone image. This kind of image is acquired with microcomputed tomography (micro-CT) to assess bone microarchitecture based chiefly on bone mineral density (BMD) measurements to improve fracture risk prediction. Disease osteoporosis can be predicted from features of CT image where a bone region may consist of several disjoint pieces. It relies on a multiresolution representation of the image by the wavelet transform to compute the multiscale morphological gradient. The coefficients of detail found at the different scales are used to determine the markers and homogeneous regions that are extracted with the watershed algorithm. The method reduces the tendency of the watershed algorithm to oversegment and results in closed homogeneous regions. The performance of the proposed segmentation scheme is presented via experimental results obtained with a broad series of images.
1. Introduction
Osteoporosis is a metabolic bone disease [1]. The disease is defined by low bone mass and microarchitectural deterioration of bone tissues leading to enhanced bone fragility. Clinically, osteoporosis is associated with an increased risk of fractures of the vertebral bone. The exact clinical estimation of bone strength and fracture risk is important for the treatment of bone diseases such as osteoporosis. It designates a major health problem in industrialized countries. This is a silent disease without any symptoms. After the age of 50, the number of osteoporotic people clearly raises [2]. A hip fracture is costly for a system of health care. Vertebral fractures are associated with a higher risk of mortality and also having another fracture in the first months following a fracture. Segmentation of bone structures from computed tomography (CT) images (shown in Figure 1) is to extract osseous tissue on the surface of bones (cortical bone, the trabecular network) from the interior region of bones (cancellous bone, space medullar).
(a)
(b)
The trabecular network detection is a very important task to help diagnose diseases such as trabecular bone loss in osteoporotic patients or in animal models of osteoporosis [1–3].
The low contrast images of the trabecular bone acquired with the sky-scan X-ray microcomputed tomography make their analysis a challenging task. Additionally, throughout acquisition, it is required to establish a high number of parameters that result in the presence of noise, nonhomogeneous illumination, and fuzzy contours. Such segmentations can be very challenging to obtain since osseous tissue does not always produce readily discernible features from soft tissue regions in CT images. Also, the trabeculae of the cancellous bone induce a certain texture pattern in the interior of the bone regions leading to large intensity variations in regions corresponding to cancellous bone.
Manual segmentation of bone regions in CT images is tedious, time consuming, and subject to observer variability; as a result, there is a strong demand for automating the segmentation procedure.
Watershed Transform [4] is a powerful morphological tool to segment images into regions of interest. This transform is able to characterize extremely composite objects. However, oversegmentation and sensitivity to noise are major problems of most watershed algorithms. The accomplishment of the Watershed Transform depends essentially on a good estimation of image gradients and the existence of obvious markers [5] for each of the objects of interest.
The automatic determination of markers is still a difficult goal to achieve. Multiscale image analysis provides a suggestion of where visually sensible objects in an image are located and also information about their relative size or importance. Wavelet transform is a good multiresolution method for edge enhancement, noise removal, and compression [6]. Morphology is also an effective method for image processing.
In this context, some methods have been proposed in the literature to analyze and simplify in different ways the image features. Although considerable research has been conducted on the segmentation problem, existing solutions to bone segmentation in CT images still remain insufficient for systematic deployment and process in a clinical setting due to drawbacks such as difficulties in algorithm initialization and parameter settings, as well as a lack of robustness to image inhomogeneity or noise. Gauch [7] has suggested an image segmentation and analysis technique employing multiscale gradient watershed hierarchies. The methodology is based on imposing a scale-based hierarchy on the watersheds associated with minima under Gaussian blurring. With the same objective, Mukhopadhyay and Chanda [8] used a hierarchy multiscale morphological operator by reconstruction for image segmentation. The result is satisfying to some extent. However, those samples of images contain a circular feature and depend on the size of the structuring element.
In another study, Kim et al. [9] proposed a multiresolution based watershed image segmentation scheme, which consists of segmenting the lowest resolution image by watershed, then using a wavelet coefficient-based region merged to reduce oversegmentation, and finally projecting the segmented low resolution image onto a full resolution image. Though, the region merging process is complicated and the calculation cost is very high.
Also, Jung [10] used a combined image denoising/enhancement technique to obtain an enhanced version of image gradient from the approximation image of the wavelet transform at the coarsest resolution based on an adaptive threshold. Then, the segmented image is projected up to higher resolutions. The method needs postprocessing to remove noise and some prior knowledge such as the rule of the postprocessing. O’Callaghan and Bull [11] proposed a combined morphological-spectral unsupervised image segmentation method. This used the subbands of the dual-tree complex wavelet transform in order to process both textured and nontextured objects in a meaningful fashion. Then, the watershed transform provides the segmented image.
Meyer [12] introduced the leveling approach, which consists of using morphological filters to remove small details of the image. The result is satisfying to some extent which has a clear edge, a small amount of noise, and a regular shape. Yet, the micro-CT trabecular bone images present a low contrast edge and intense amount of noise.
In this paper, we propose to segment the trabecular bone images using a watershed method with a multiresolution analysis to compute a gradient recognition and controlled markers. In the first step we suggest creating a reconstructed morphological gradient image by using the multiresolution representation of the image which reduces the oversegmentation of the watershed algorithm. Then, we proposed using the wavelets transform to locate the regions of interest and the morphological operations to obtain and preprocess the skeleton in order to obtain markers that allow robust segmentations in images.
The paper is organized as follows: in Section 2, the watershed and the wavelet transform are described; Section 3 focuses on our approach for trabecular bone segmentation and its implementation. Before concluding, the results are discussed in Section 4.
2. The Watershed Transform and the Wavelet Transform
2.1. The Watershed Transform
The watershed transformation is a morphological tool dedicated to image segmentation [4, 5, 13]. According to a flooding paradigm, this transformation considers each gray level (or color or texture) as topographic data as shown in Figure 2.
The principle is to first search out all the local minima where to perforate holes. Then, immerse the model into a lake at a uniform speed as the time. During the flooding stage, the water level slowly rises, covering an increasing amount of the image. To prevent the fusion of two catchment basins corresponding to different regional minima, a dam is erected between them at the points where they congregate each other. The watershed line corresponds to a series of this constructed dam and local extremes. For digital images, watershed constitutes closed contour of one-pix width [6]. Different approaches are employed to use the watershed principle for image segmentation. One of the most common watershed algorithms was introduced by Meyer [6]. This flooding process is performed on the gradient image. The implementation of this algorithm uses the hierarchical queue [5].
This transformation is characterized by its fast computing and high accuracy in locating the weak edges of adjacent regions. But normally this will lead to an oversegmentation of the image. Therefore, many researchers [14–18] have proposed various methods for improvement, such as a preprocessing step before watershed transformation.
2.2. The Wavelet Transform
The wavelet transform [19] is a mathematical tool that can represent signal (image) in both spatial and frequency domains. Therefore, it is used to represent images in multiple resolutions. The wavelet transform decomposes signals into basic functions that are dilations and translations of a single prototype wavelet function: where are obtained by translating and dilating the wavelet function . The discrete wavelet transform coefficient is the estimation of signal components in the time-frequency plane. This coefficient is obtained by this product:
According to Mallat’s pyramid algorithm [20], the signal is decomposed using a low-pass filter “” and a high-pass filter “” and each filter bank is then sampled at a half rate (1/2 downsampling) of the previous frequency. By repeating this procedure, it is possible to obtain wavelet transform of any order. In the case of an image, the filtering is implemented in a separable way by filtering the lines and columns. According to this algorithm, the original image can be transformed into four subimages. A comprehensive survey of the wavelet transform is provided in [21]. This pyramidal representation is not desirable for estimation or detection problems. In order to preserve the invariability by translation, the downsampling operation must be suppressed and the decomposition obtained is called a stationary wavelet transform [21]. In practice, the structure in cascade of the filter bank does not change; the only operation of downsampling is suppressed. In this work, we choose the Haar wavelet [19, 20], because it is orthogonal, has small support (for this reason, it provides good localization in space), and requires small computational complexity (linear with respect to the size of the input image).
3. The Proposed Method
The proposed algorithm can be summarized into the following steps. (i)Select a desired scale , and compute the stationary wavelet transform of the image up to this scale. (ii)Create a reconstructed gradient image by combining the morphological gradient magnitude found at each scale, as shown in Figure 3. (iii)Localize the foreground objects using stationary wavelet reconstruction of the only high frequency part and threshold the obtained image. (iv)Clean the edges of this binary image using morphological reconstruction opening and closing operators. (v)Calculate the skeleton of this preprocessed image to obtain the perfect suitable markers. (vi)Compute the background markers. (vii)Superimpose these markers in the reconstructed gradient image. (viii)The latest obtained image is given as input to the watershed transform.
These steps are detailed next.
3.1. Multiresolution Morphological Gradient
The gradient magnitude [22] is used often to preprocess a gray-scale image prior to using watershed transform. The combination of the morphological gradients in different resolutions [23] is insensitive not only to noise but also to extract various fineness of the edges, because from high-resolution images, small details can be detected, while from low-resolution images, larger structures are detected. Thus, to obtain a reconstructed gradient image, the original image is transformed into a subimage using stationary wavelet transform. We retained only the low frequency part of this subimage. The morphological gradient magnitude is calculated for each resolution as where and stand, respectively, for dilation and erosion. “” is the structuring element ( () rectangular SE). The reconstructed gradient is obtained by combining the gradient magnitude images found at each scale such as is the coarsest resolution.
3.2. Creating Markers
The oversegmentation of the watershed algorithm can be reduced by finding markers for each object. However, automatic marker selection is a difficult process and in most applications human interaction is required. In addition, if inappropriate markers are selected, the final result changes dramatically. Based on the analysis of the characteristics of the trabecular bone image, the traditional function that detects the regional minimum as markers is not efficient in the trabecular bone. This method identifies many minimums in the same region, so we obtain an over-segmentation in this region after the segmentation. Therefore, a novel combination method is proposed to detect markers. Firstly, wavelet transform and morphology are used to locate the elements and get the subimages including the elements. The major processes of this step are described as follows. (i)Get the input image. (ii)Adopt “Haar” wavelet for the wavelet transform of the input image. (iii)Set the low frequency component after wavelet decomposition as 0; the high frequency part does not change. (iv)Wavelet reconstruction: get the reconstructed image. (v)Threshold the reconstructed image using Otsu’s method [14] and binary image is obtained.
Figure 4 shows the localization of the object of interest.
(a)
(b)
Based on the investigation of the characteristics of the micro-CT trabecular bone image, we suggest computing the skeleton of the objects of interest as markers, provided that each object only contains one marker. In order to create foreground and background markers, a binary image enhancement is needed. So, firstly, we use the morphological opening by reconstruction followed by the morphological closing reconstruction to remove small particles due to noise. These operators [6] are listed as follows: Let be the two images defined on the same domain and . (circular SE) is the flat structuring element.(i)Opening by reconstruction: the reconstruction of from is obtained by iterating elementary geodesic erosions of above until stability is reached and denoted as where can be obtained by iterating “” elementary geodesic erosion and the geodesic erosion is defined as (ii)Closing by reconstruction: the reconstruction of “” from “” is obtained by iterating elementary geodesic dilations of “” under “” until stability is reached and denoted as where can be obtained by iterating “” elementary geodesic dilation and the geodesic dilation is defined as
Then, to get smooth edge foreground objects of the binary image, we use some morphology operations such as closing and dilating. These morphology operations fill the gaps between pixels and smooth their outer edges. The size of the gaps between pixels must be the maximum of structuring element size.
Figure 5 shows the result of binary image enhancement with morphology filters.
(a)
(b)
The foreground and background markers are created by calculating the morphology skeleton of the preprocessed binary image. In order to find the skeleton of trabecular network or medullar space in trabecular bone image, we use the algorithm of Gonzales and Woods [24] for thinning binary regions. Region points are assumed to have value 1, and background points value 0. Thinning is a morphological operation, somewhat like erosion or opening. It has-been designed as an iterative edge deleting items on the area. By considering the constraints that deletion of these points [24, 25] does not remove end points nor break connectivity nor cause excessive erosion of the region.
The thinning [24] of an image by a structuring element is related to the hit-and-miss transform and it is defined as where the subtraction is a logical subtraction defined by
The hit-and-miss transform of an image by structuring elements (circular SE, size 2) and (circular SE, size 5) (they should not be having any overlapping element) is defined as where and stand, respectively, for erosion and the complement of “.”
Figure 6 shows the obtained foreground and background markers.
(a)
(b)
3.3. Watershed Segmentation
Before we apply the watershed transformation to the reconstructed morphological gradient, we impose foreground and background markers as minimums. Then, the flood will come off only from these marker pixels. The modified gradient with imposed minima is given by
The result image is given as an input to the watershed algorithm. Let us represent the watershed segmentation operator by ; the output image is defined as
4. Experimental Results
The proposed segmentation algorithm is applied to a wide number of micro-CT trabecular bone images of human bone biopsy, taken from the femoral neck [2]. The results are subjectively analyzed. The trabecular network is effectively segmented from the medullar space. The markers are selected adaptively, depending on the characteristics of the trabecular bone image. They are employed to join adjacent regions because they should be corresponding to the same object and should eliminate the split of the wanted objects. Although the proposed algorithm prevents oversegmentation, it preserves the subfeatures of the objects required for getting morphological calculation results and reconstructed 2D slices. These segmentation results can be used in the feature extraction and classification stages to distinguish a normal bone from an osteoporotic bone. The segmentation results as produced by the proposed algorithm are shown in Figure 7.
(a)
(b)
(c)
4.1. Comparison with Other Techniques
To evaluate the performance of the presented method, we compare the proposed technique with other preprocessing techniques known to improve the watershed segmentation, namely, Marker-Controlled Watershed Segmentation [24], combined morphological-spectral unsupervised image segmentation method (VPSN) [11], and watersheds segmentation method using wavelets to enhance the image [16–22]. Results were analyzed in both quantitative metric and visual inspection.
The segmentation results using the classical Marker-Controlled Watershed, the method combining the watershed with the dual-tree complex wavelet transform (VPSN), and the method using the wavelet to enhance the image after applying the watershed are shown in Figures 8(c), 8(d) and 8(e). Comparing with Figure 8(b), it is obvious that our proposed strategy produces a much more accurate segmented image. Also, we applied our strategy to other low contrast images such as Figures 9 and 10.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
Figures 9(b) to 9(e) and Figures 10(b) to 10(e) show the segmentation results obtained by the proposed technique and the other proposed techniques known to improve the watershed segmentation cited above. In fact, visual inspection shows that the proposed technique produced the best result.
Figure 11 compares computational time (in seconds) of the four segmentation methods of the industrial circuit image.
The proposed technique using the multiresolution to estimate the reconstructed morphological gradient image and the high frequency subbands to locate and compute the markers is significantly faster than the other methods.
Furthermore, results were evaluated in quantitative points of view. The chosen quantitative metric was a percentage error obtained by computing the difference between the segmented image and the reference image according to the equation below: where “” represents a binary image of the reference image, “” represents a binary image of the segmented image, and the “” represents the number of the pixels of the corresponding binary image.
For example, a quantitative analysis (according to (14)) indicates that segmentation errors (in %) corresponding to Figures 8(b) to 8(e) are, respectively, 12.73, 17.56, 65.74, and 22.56% (reference image was obtained by manual delineation). The proposed technique produced the smallest segmentation error. Also the segmentation errors (in %) corresponding to Figures 9(b) to 9(e) are, respectively, 20.47, 27.50, 25.22, and 84.27%. Once more, the proposed technique produced the smallest error.
5. Conclusions
A novel method using morphological operations, wavelet, and watershed segmentation is proposed to automatically detect the trabecular network from micro-CT trabecular bone images. The proposed method does not require human intervention. The experiment shows that the proposed algorithm produces accurate boundaries of the region of interest. In addition to that, the algorithm is implemented automatically leading to efficient computation work. We should also notice that the segmentation result depends on the quality of the gradient image and the suitable markers.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work has been conducted in the Research Unit Sciences and Technologies of Image and Telecommunications (SETIT) ISBS.