Abstract
Automatic segmentation and three-dimensional reconstruction of the liver is important for liver disease diagnosis and surgical treatment. However, the shape of the imaged 2D liver in each CT image changes dramatically across the slices. In all slices, the imaged 2D liver is connected with other organs, and the connected organs also vary across the slices. In many slices, the intensities of the connected organs are the same with that of the liver. All these facts make automatic segmentation of the liver in the CT image an extremely difficult task. In this paper, we propose a heuristic approach to segment the liver automatically based on multiple thresholds. The thresholds are computed based on the slope difference distribution that has been proposed and verified in the previous research. Different organs in the CT image are segmented with the automatically computed thresholds, respectively. Then, different segmentation results are combined to delineate the boundary of the liver robustly. After the boundaries of the 2D liver in all the slices are identified, they are combined to form the 3D shape of the liver with a global energy minimization function. Experimental results verified the effectiveness of all the proposed image processing algorithms in automatic and robust segmentation of the liver in CT images.
1. Introduction
Liver diseases have become one of the most common causes of deaths in the world. Researchers have focus on the prevention and treatment of liver diseases for many years. In recent years, computed tomography (CT) imaging has been widely used in liver disease diagnosis and surgical treatment because tumors or hepatic lesions could be observed easily from the CT image. For the captured CT images, the liver slices need to be examined in the two dimensions one by one. Consequently, it lacks an overall image of the 3D liver. Furthermore, it takes clinicians considerable time to view all the slices and diagnose the disease or evaluate the liver function based on the information divided and presented in different images. Therefore, it is desirable that 2D liver slice is segmented from CT images, and the 3D liver is reconstructed automatically and robustly beforehand. Thus, the clinicians could get the information of the 3D liver at a glance and diagnose the liver disease or evaluate the liver function more conveniently. This desire has led researchers worldwide to devote themselves to the research of coming up with automatic and robust liver segmentation methods. After so many years of research, it remains an open problem because the liver is adjacent to many other organs, such as the kidneys, spleen, stomach, intestines, and bones. In many cases, the intensity of the liver and that of the adjacent organ is indistinguishable. In addition, the shape of the liver varies according to the individuals. As a result, automatic and robust segmentation of the liver from the CT images remains as one of the most challenging artificial intelligence task for many decades.
After so many years of study and research, quite a few liver segmentation methods have been proposed, though none of them have achieved adequate accuracy so far. Among them, methods based on statistical and probabilistic models became the most popular ones [1–4]. Yet, such kinds of methods require a large size of training samples, which decreases the segmentation efficiency tremendously. Even if dictionary and sparse coding is used to reduce the redundant information, the segmentation efficiency is still not satisfactory. In recent years, deep learning and convolutional neural network have been used to segment the liver in the CT images, and the reported accuracy appears to be promising [5, 6]. Similarly, deep learning or convolutional neural network-based methods rely heavily on the training datasets to yield accurate segmentation results. In other words, the accuracy will not be acceptable if the training datasets are not similar enough to the tested case. In reality, the liver of the patients varies tremendously, and it is thus difficult to acquire a complete training datasets. In such situations, both statistical/probabilistic models-based methods and neural network-based methods might fail completely in identifying the livers of some individuals.
There are also many other methods that do not rely on training datasets. The popular segmentation methods used in identifying the boundary of the liver in CT images include active contour [7], threshold selection [8], level set [9], and graph cuts [10–12]. Since its surrounding tissues and connected organs in different slices are different, one segmentation method or several combined segmentation methods [7–12] usually could not segment the liver robustly and fully automatically. Hence, heuristic methods are proposed to segment the liver, its surrounding tissues, and its connected organs as a whole [13–15]. For instance, the surrounding tissues are segmented first and then serve as a constraint for the segmentation of the liver. In this way, the surrounding tissues are changed from interference factors to the markers of constraint. Thus, the liver could be segmented robustly and automatically [16].
In this paper, we also propose a heuristic method based on multiple thresholds selection and morphological operations. With a global threshold, we could robustly segment all the connected organs in each CT slice, such as the ribs, the spine, the heart, the kidney, and the stomach. With another global threshold, we could robustly segment the liver with the logic not operation of the segmentation result of the connected organs. Therefore, robust threshold selection becomes a critical step in the proposed method. We utilize the slope difference distribution-based threshold selection (SDDTS) method to calculate multiple thresholds in this research work. The robustness of the slope difference distribution-based threshold selection method and its advantage over state-of-the-art segmentation methods have been verified in our previous studies [17, 18]. At first, we calibrate the parameters of the slope difference distribution-based threshold selection method with several typical CT images and then use the calibrated parameters for all the CT slices. After segmentation by the global threshold, the segmentation result is filtered by minimizing the Gibbs energy function [19] to reduce inhomogeneity. Then, morphological operations [20] are used to merge divided parts of the same organ, and the merged organ is used as liver segmentation constraint. Since the intensities of the surrounding tissues and that of the liver might be completely the same, we use centroids of the segmented ribs and the spine to fit a curved line, which is then used to separate the liver and the surrounding tissues. After liver segmentation, its boundary is further refined by spline filtering.
2. Methods
The proposed method is heuristic, and it contains a series of image processing algorithms that vary depending on the index of the CT slice. The core image processing algorithms include (1) threshold selection based on slope difference distribution [17, 18] and image segmentation with the selected threshold; (2) energy minimization of segmentation result to eliminate noise; (3) morphological filtering; (4) morphological merging; and (5) spline filtering. The flowchart of the proposed method is shown in Figure 1. The multiple thresholds are calculated from the inputted 2D CT slice. One threshold is used for liver segmentation, and the other thresholds are used for constraint segmentation. The core image processing algorithms are used during constraint segmentation and computing the boundary of the segmented liver.
2.1. Threshold Selection
The slope difference distribution is computed from the histogram of the image, and it reflects the global variation rate of the histogram. With the assumption that the thresholding point between two pixel classes varies greatest, the thresholding point could be computed based on the variation rate of the histogram, i.e., the slope difference distribution. The slope difference distribution is formed by a series of slope difference that is calculated at each sampled point of the smoothed histogram. At each sampled point of the smoothed histogram, points on its left are selected to fit a line model and points on its right are selected to fit another line model. The slopes of these two line models are calculated. The left slope is subtracted from the right slope, and the slope difference at this sampled point is obtained. First of all, the histogram of the image needs to be calculated and smoothed as follows.
The grey-scale values of the image are scaled to the interval [1, 255], and then the histogram distribution is calculated as follows:where denotes the total number of pixels that equals to and denotes the maximum number of pixels that occurs at in the interval [1, 255]. The histogram distribution, , is then transformed to the frequency domain by the discrete Fourier transform (DFT):
Only the low frequencies from 0 to and the symmetric frequencies from to 255 are kept:where is the bandwidth of the low-pass DFT filter. Its value and the value of are obtained by parameter calibration [17, 18]. In general, there are always optimal values of and to yield the optimal threshold. However, it is still challenging to determine the optimal and by simple parameter calibration. When the values of and are not determined properly, the threshold will be determined inaccurately. The reason we choose the Fourier transform-based filter instead of other popular filters is based on the quantitative evaluation and comparison [19]. We found that the Fourier transformation-based filter outperforms other filters significantly in this conducted research work. Transforming from the frequency domain back to the spatial domainwhere is the smoothed histogram distribution.
To compute the slope difference, we fit two-line models on both sides of the sampled point. The line model is formulated as follows:where , is the sampled point on the smoothed histogram distribution and is the slope of the line model. The coefficient of the line model is computed as follows:where is the design matrix of the least square fitting method and is the input data vector. Moreover, are the adjacent points at the left side of the point , and are the adjacent points at the right side of the point . The left slope and the right slope, and at point , are then obtained from Equation (6). The slope difference at point is then computed as follows:
The discrete version is denoted as , and its continuous version is denoted as that is named as slope difference distribution. Let the derivative of equal zero and solve it, the valleys with greatest local variations are obtained. The positions where these valleys lie are the thresholds that separate different pixel classes. One fundamental property of the slope difference distribution is that the positions of the valleys change monotonically with the line model fitting parameter . Hence, the parameter could be calibrated to yield the optimal threshold. After the optimum threshold is selected, the image is segmented by the following equation:where is the original image and is the binarized image. is the index of the pixel in the image. Since multiple thresholds are needed to segment different organs, Equation (9) is the basic format of segmenting a single object. More conditions need to be added in Equation (9), when several independent segmentation results are combined to segment the liver.
2.2. Energy Minimization
There are many popular noise reduction methods, e.g., the discrete Fourier transformation-based filter [20] and the wavelet image processing [21]. The noise reduction methods are usually applied beforehand in the preprocessing [20]. Since we did not reduce the noise beforehand, we apply a noise reduction procedure immediately after the segmentation. The noise would cause region inhomogeneity, and this inhomogeneity could be formulated by the Gibbs distribution:where is the potential function associated with clique . The clique is defined as a set of sites such that any two elements in the clique are neighbours of each other [19]. is the total number of pixel classes. is defined by the following equation:where is the constant and 1 is its default value. To reduce the noise in the binarized image, the total energy over the whole image is minimized as follows:where is the filtered image with Gibbs energy minimization; and are the pixel indexes of the image in the vertical direction and in the horizontal direction, respectively; and and are the resolution of the image in the vertical direction and in the horizontal direction, respectively.
2.3. Morphological Filtering
Besides the noise, there are many small binarized blobs that do not belong to the liver or the segmented organ. These small blobs could not be removed by energy minimization. A popular way to remove these small blobs is to count the areas of all the blobs and remove some of the blobs based on their areas morphologically. However, there are also many situations in which these small blobs are connected with the liver or the segmented organ, which make the removal of them more difficult. To remove all these interference blobs, we propose a morphological filtering method that contains the following steps.
Step 1. Erode the segmentation result morphologically as follows:where is the 4-connected structure element with the disk shape and its radius is 1, is the point in the structuring element , and is the translation vector.
Step 2. Repeat Step 1 times. The default value of is 8.
Step 3. Dilate the segmentation result morphologically as follows:where denotes the symmetric or supplement of .
Step 4. Repeat Step 3 times.
In summary, the proposed morphological filtering method removes the small blobs by a repeating morphological erosion process first. Then, it restores the eroded liver or other organs by a morphological dilation process with the same repeating times.
2.4. Morphological Merging
On the contrary, there are situations where the segmentation result of the organ specially the stomach is split into different parts. To utilize the segmentation result more effectively, it is required that the split parts are merged into a united one. To unite these split parts, we propose a morphological merging method that contains the following steps.
Step 1. Dilate the segmentation result morphologically as follows:
Step 2. Repeat Step 1 times. The default value of is 16.
Step 3. Erode the segmentation result morphologically as follows:
Step 4. Repeat Step 3 times.
As can be seen, the morphological merging is the opposite process of the morphological filtering. It connects the split parts and merges them into one united one by a repeating morphological dilation process first. Then, it restores the dilated organ by a morphological erosion process with the same repeating times.
2.5. Spline Filtering
After the liver is segmented, its boundary is extracted first. The final smooth boundary is computed by minimizing the energy between the extracted boundary and the fitted polynomial spline by the following equation:where is the smoothing factor and its default value is 0.5.
2.6. Liver Segmentation
We use a typical CT image as shown in Figure 2(a) to illustrate the proposed segmentation method. First, the thresholds are computed with the slope difference distribution as shown in Figure 2(b). Four thresholds are calculated, and they are shown in Figure 2(b), respectively. The threshold to segment the spine and ribs is denoted by the green asterisk in the green circle, the threshold to segment the body is denoted by the blue asterisk in the green circle, the threshold to segment the stomach is denoted by the red asterisk in the green circle, and the threshold to segment the liver is denoted by the black asterisk in the green circle. With the threshold , the spine and the ribs are segmented as follows:where denotes the segmented spine and ribs, and it is shown in Figure 2(c). The centroids of these blobs are calculated as the means of all the pixels they contain. Based on the computed centroids, the blobs are divided into two classes, the first class of blobs that include the blobs on the left and the second class that includes the blobs on the right and around the center. The blobs in the second class are filtered based on their areas, and only the largest one (the spine) is kept. The centroids of the blobs in the first class and the largest one in the second class are used to fit a second-order curve by the least squared method (Equations (6–8)) as shown in Figure 2(d). This fitted curve is used to separate the tissues with almost the same intensity from the liver. Figure 2(e) shows the fitted curve overlaying on the original image . As can be seen, the liver and its neighboring tissue on the bottom are separated successfully.
With the threshold , the whole body is segmented as follows:where denotes the segmented body, and it is shown in Figure 2(f). Figure 2(g) shows the segmentation result of the body after morphological filtering. As can be seen, all the small blobs are removed successfully. The filtered body is then eroded morphologically (Equations (14)–(16)) by the disk structuring element with the radius of . is calculated as the average width of the segmented ribs, and it is 16 in this specific example. The eroded body is then subtracted from the filtered part, and the constraint circular part is obtained as shown in Figure 2(h). This constraint circular part is also used to eliminate the surrounding tissues that have the same or similar intensities with that of the liver. The constraint circular part overlaying on the original image is shown in Figure 2(i). As can be seen, it is adjacent to the liver, but almost not covering any liver part.
With the threshold , the stomach is segmented as follows:where denotes the segmented stomach, and it is shown in Figure 2(j). It is then filtered by energy minimization, and the result is shown in Figure 2(k). The filtered result is subtracted by the segmented bones shown in Figure 2(c), and the subtraction result is shown in Figure 2(l). The morphological filtering is applied on the subtraction result, and then, only the parts on the right of the spine are kept as shown in Figure 2(m). The morphological merging is applied on the kept parts, and the merged stomach is shown in Figure 2(n). With the constraint of fitted curve , the circular part , and the merged stomach , the liver is segmented by the threshold as follows:where denotes the segmented liver, and it is shown in Figure 2(o). It is then filtered by energy minimization, and the filtered result is shown in Figure 3(a). As can be seen, there are still significant interference blobs. Immediately following the energy minimization, the morphological filtering is applied, and the filtering result is shown in Figure 3(b). As can be seen, the interference blobs are reduced significantly, and the left blobs could be removed by the morphological area filtering. The largest blob (the liver) is kept after the morphological area filtering and it is shown in Figure 3(c). The boundary of the segmented liver is extracted and smoothed by spline filtering. Figure 3(d) shows the smoothed liver boundary overlaying on the original image. As can be seen, the extracted boundary matches the liver very well.
(a)
(b)
(c)
(d)
2.7. Three-Dimensional Reconstruction
After the two-dimensional livers in all the CT slices are segmented and their boundaries are extracted, each boundary is sampled evenly with the same number of points; i.e., the distance between any two adjacent sampled points on the same boundary is almost the same. The number of the sampling points is chosen as 200 in this study. Then, all the sampled points from all the slices are stacked together according to their practical pixel distances during the CT scanning. As a result, the three-dimensional coordinates of all the sampled points are obtained. The index of the sampled points in each slice is aligned with each other in the stacking direction (z direction). Then, we get 200 curved lines in the stacking direction, and the number of the points on the curved line equals the stacking number . We resample the curved lines by the following spline interpolation filter:where denotes the smoothed point and denotes the original sampled point. is a smoothing factor and is the fitted spline function ().
3. Results and Discussion
Since the shapes of the liver and the surrounding organs also vary significantly in different slices, the number of the pixel classes in different slices varies accordingly. The organ that is adjacent to the liver is not fixed, and it might be the stomach, the kidney, or the heart. We thus define the liver segmentation into three cases depending on its adjacent organ. In the first case, there is only stomach adjacent to the liver. In the second case, there is kidney adjacent to the liver. In the third case, there is heart adjacent to the liver. In different cases, the number of the pixel classes is adjusted automatically based on the detected thresholds by the slope difference distribution. Then, the five core image processing algorithms are applied one by one. The optimal values of the input parameters are determined for different cases, respectively. After the optimal parameters are calibrated by trial and analysis, they will be used for segmenting the livers in similar CT datasets. For the CT datasets with significant differences, the input parameters should be recalibrated again. As a result, it might require significant manual intervention to determine the optimal parameters for a specific-type CT dataset before the proposed method is run fully automatically for the whole CT dataset. The average time to process one CT image is 2.39 seconds in MATLAB with the i7-6700 CPU. Some typical liver segmentation results for the first case are shown in Figure 4. Some typical liver segmentation results for the second case are shown in Figure 5. Some typical liver segmentation results for the third case are shown in Figure 6. As can be seen, all the segmentation results are acceptable for clinical usage. With all the segmented two-dimensional liver from different slices, the whole three-dimensional liver is reconstructed by the method described in the above section. Figure 7 shows the reconstructed three-dimensional liver.
(a)
(b)
We use two public datasets to evaluate the proposed method further. The first public dataset (https://eee.deu.edu.tr/moodle/mod/page/view.php?id=7872) is liver transplantation donor database that provided by Emre Kavur from Dokuz Eylul University with the permission of Dokuz Eylul University Hospital. In this set, there are 20 upper abdominal CT image series that belongs to different patients and six sets are training sets. Only ground truths of the training sets are provided. Since the proposed method does not require training, we use the provided training sets to evaluate our method. The computed mean ratios of the volumetric overlap error (VOE), relative volume difference (RVD), average symmetric surface distance (ASD), root mean square symmetric surface distance (RMSD), and maximum symmetric surface distance (MSSD) were 9.6 ± 2.2%, 4.2 ± 2.5%, 1.7 ± 0.9 mm, 2.4 ± 1.1 mm, and 9.2 ± 3.1 mm, respectively. We also show some qualitative results in Figure 8. As can be seen, the automatic identified boundary of the liver matches the manually identified boundary well.
The second public dataset, 3D-IRCADb01 (https://www.ircad.fr/research/3d-ircadb-01/), consists of 20 CT scans with corresponding ground truth provided by IRCAD, the French Research Institute against Digestive Cancer. We show the quantitative comparisons with state-of-the-art methods [22–25] in Table 1. As can be seen, the comparisons are favorable. The proposed method does not require any work for training as state-of-the-art methods [22–25] do while the achieved accuracy is similar or better.
The accuracy of the proposed method differs significantly across different CT liver datasets because the proposed method is heavily relying on the accuracy of threshold selection. Yet, the optimal accuracy of slope difference-based threshold selection could not be guaranteed currently during unsupervised segmentation. In the near future, we will research on algorithms to make the slope difference-based threshold selection method be capable of always selecting optimal thresholds without supervision.
4. Conclusions
In this paper, a heuristic approach is proposed to segment the liver in CT images fully automatically. It calculates multiple thresholds simultaneously based on the slope difference distribution and then segment the CT image into different meaningful regions with these automatically calculated thresholds. The liver is segmented robustly with the constraint of the surrounding organs or tissues segmented beforehand. To reduce the noise, Gibbs distribution is utilized to minimize the global energy. A morphological filter is proposed to remove the small interference blobs. A morphological merging method is also proposed to unite the divided parts. Experimental results verified all the proposed image processing algorithms and the proposed approach. Since the proposed approach is very efficient and does not require any training datasets, it is thus very promising for clinical usage.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
The authors thank the hospital at Qingdao University for providing the dataset.