Patient-specific probabilistic atlas combining modified distance regularized level set for automatic liver segmentation in CT

Abstract

Liver segmentation from CT is regarded as a prerequisite for computer-assisted clinical applications. However, automatic liver segmentation technology still faces challenges due to the variable shapes and low contrast. In this paper, a patient-specific probabilistic atlas (PA)-based method combing modified distance regularized level set for liver segmentation is proposed. Firstly, the similarities between training atlases and testing patient image are calculated, resulting in a series of weighted atlas, which are used to generate the patient-specific PA. Then, a most likely liver region (MLLR) can be determined based on the patient-specific PA. Finally, the refinement is performed by the modified distance regularized level set model, which takes advantage of both edge and region information as balloon force. We evaluated our proposed scheme based on 35 public datasets, and experimental result shows that the proposed method can be deployed for robust and precise liver segmentation, to replace the tedious and time-consuming manual method.

Keywords:

• Probabilistic atlas
• patient-specific
• level set
• liver segmentation

1. Introduction

Abdominal liver segmentation is a fundamental step in liver quantitative analysis. However, due to the low contrast and variable shape, robust and accurate liver segmentation is still challenging, and it is often insufficient to segment the liver from CT using traditional segmentation techniques. Currently, varieties of approaches have been proposed for liver segmentation [1]. These studies segment liver either based on pure image information or based on model.

For pure image information-based approaches, Ruskó et al. [2] first used the 3 D region growing to initialize the liver contour, and then utilized post-processing techniques, such as separating the heart, removing the inferior vena cava to refine the result. However, the grayscale-based liver segmentation results are often not satisfactory, since the liver and the other adjacent organs have similar gray values in CT.

In the case of clustering or classification-based approaches, novel fuzzy clustering methods proposed by Jiang et al. (e.g. Co-FCM [3] and KL-TFCM [4]), have significantly improved the performance of segmentation. Li et al. [5] proposed to integrate the both probabilistic shapes and global intensities into the Bayes classification frame, to achieve simultaneous segmentation of multiple organs from the abdominal CT. Meanwhile, multi-feature classifier based multi-layer liver segmentation method is explored by Selver et al. [6].

Liver segmentation is also performed by the snake model. Lu et al. [7] proposed a semi-automatic snake-based method. The method first initializes the snake model by the existing segmentation result (3 D mesh), and in the subsequent deformation process, the internal force, the external force and the newly introduced attractive force are used together to determine the action of each node on the grid surface.

Compared with the snake model-based approaches, the level se is another type of deformation model-based method, which is more commonly employed for liver segmentation. Generally, the level set-based methods tend to combine with other information (e.g. statistical prior knowledge, local or global image features) to improve the accuracy and reduce the computational complexity. Yang et al. [8] proposed a semi-automatic segmentation method for liver CT, which uses level set to achieve refinement. Kohlberger et al. [9] presented a level set approach utilizing local shape and appearance prior model, in which the local intensity statistics and local curvature information are integrated, while tracking the evolution of each level set point, and thus achieves a satisfactory result.

Since the gray scales of the organs in the abdominal CT are similar, integrate the shape prior information may significantly improve the precision of the segmentation. Zhang et al. [10] utilized the statistical shape model combining with the optimal surface acquisition to achieve the automatic liver segmentation. Tomoshige et al. [11] proposed a method based on conditional statistical shape model. The method first uses the maximum a posteriori (MAP) estimation to extract conditional features, and then by employing the level set to obtain a liver shape prior estimation. Finally, an accurate refinement of the liver contour was achieved by employing graph cut method.

Atlas-based method is another popular approach for robust liver segmentation. Zhou et al. [12] first used liver probabilistic atlas for automatic liver segmentation in CT. The probabilistic atlas constructed in this method provides two kinds of liver information: liver position and liver grayscale, and experiment shows that their atlas-based method can reduce the calculation cost of liver segmentation process, and can significantly improve the robustness.

Although the accuracy of liver segmentation is greatly improved by the presented methods, the existing methods are still difficult to meet the requirements in clinical practice. Aiming at the low contrast between organs in the abdominal CT image, organ pathology and the difference of organ shape between individuals, this paper proposes an automatic liver segmentation method based on the combination of probabilistic atlas and modified distance regularized level set model.

2. Method

In this section, the proposed segmentation framework will be described in detail. The flowchart is depicted in Figure 1.

Figure 1. The flowchart of the framework.

2.1. Denoising

For noise reduction, this paper employs the Catte_PM anisotropic diffusion filter [13] to denoise the liver CT image. Figure 2 shows a result of filtering the liver CT image using the Catte_PM model. As seen in the figure, the Catte_PM model can not only smooth the inner region of the organ, but also effectively preserve the edge information between organs.

Figure 2. Denoising of liver CT image.

2.2. Patient-Specific probabilistic atlas

For a given test dataset (patient image), the Patient-Specific PA can be dynamically generated by the following steps. First, the image similarity of all the maps and the test set is calculated, and the atlases are arranged in descending order in accordance with the obtained image similarity. The atlas is composed of N training CT datasets with corresponding “golden standards”. The similarity between the test patient and training atlas can be obtained according to the Normalized Cross Correlation (NCC) [14] method, and by employing the method proposed by Chu [15], a weighted probabilistic atlas is calculated.

2.3. Most likely liver region (MLLR)

After constructing the weighted PA of the testing image, the PA is then thresholded, and a binarized liver parenchymal region of interest (ROI) is generated, which is called the most likely liver region (MLLR), as shown in Figure 3. In order to minimize the estimation error of the MLLR, a morphological dilation operation is performed after binarization, in which spherical structural elements are used. The resulting segmentation is then assigned as the initial contour for the following level set refinement.

Figure 3. The Most likely liver region (MLLR).

2.4. Modified distance regularized level set

We descript a Modified Distance Regularized Level Set-based refinement method, in which the evolution is driven by balloon forces of edge and region information, and thereby the active contour can evolve to the edge of the liver with sufficient balloon force or sufficient number of iterations. Moreover, edge leaks do not easily occur, resulting in more accurate results.

The basic edge-based Distance Regularized Level Set model utilized in this paper is proposed by Li [16]. Given level set equation φ: image domain Ω→R, and an energy function ε(φ) can be defined as follows (1) $$\epsilon \left(\mathit{\phi }\right)=\mu D_p\left(\phi \right)+{\epsilon }_{\mathit{ext}}\left(\phi \right)$$(1) where μ > 0. D_p(φ) is a penalty term, which is aimed to avoid the re-initialization of the level set, and it is defined as follows: (2) $$D_p\left(\phi \right)={\int }_{\Omega }p\left(\left|\nabla \phi \right|\right)dx$$(2) (3) $$p\left(\left|\nabla \phi \right|\right)=\{\begin{array}{l}\frac{1}{{\left(2\pi \right)}^2}\left(1-\hspace{0.17em}\text{cos}\hspace{0.17em}\left(2\pi \left|\nabla \phi \right|\right)\right),\nabla \phi \le 1\\ \frac{1}{2}{\left(\left|\nabla \phi \right|-1\right)}^2,\nabla \phi \ge 1\end{array}$$(3) where p(|▽φ|) is a double-well potential function. The main role of this function is to allow the binary Heaviside function to be directly initialized for the level set, and the function ε_ext(φ) is an external force that can be defined according to the edge or the regional information to drive the evolution of the level set.

An edge-based DRLSE level set was constructed by adding the edge information ε_ext(φ) of image. Define the edge indication function g: (4) $$g=\frac{1}{1+{\left|\nabla G_{\sigma }\ast I\right|}^2}$$(4) where G_σ represents the Gaussian kernel function with variance σ, and ε_ext(φ) is defined as follows： (5) $${\epsilon }_{\mathit{ext}}\left(\phi \right)=\lambda L_g\left(\phi \right)+\alpha A_g\left(\phi \right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\hspace{1em}=\lambda {\int }_{\Omega }g\delta \left(\phi \right)\left|\nabla \phi \right|dx+\alpha {\int }_{\Omega }\mathit{gH}\left(-\phi \right)dx$$(5) where Ω is the image domain, λ andαare both constant coefficients, δis the Dirac function and H represents the Heaviside function, L_g(φ) represents the weighted length of the curve and A_g(φ) is the evolution of the balloon force and control curve.

In our modified level set model, a new balloon force based on regional information （Signed Pressure Force, SPF）proposed by Zhang [17] is combined into edge-based DRLSE, and SPF(I(x)) is defined as follows, (6) $$\mathit{SPF}\left(I\left(x\right)\right)=\frac{I\left(x\right)-\frac{c_1+c_2}{2}}{\max \left(\left|I\left(x\right)-\frac{c_1+c_2}{2}\right|\right)}$$(6) in which c_1, c₂ are the average gray-level constants inside and outside the curve, which can be obtained by solving the C-V model [18]. As the SPF(I(x)) is obtained, it can be incorporated into the balloon force of Equation (5)(5) $${\epsilon }_{\mathit{ext}}\left(\phi \right)=\lambda L_g\left(\phi \right)+\alpha A_g\left(\phi \right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\hspace{1em}=\lambda {\int }_{\Omega }g\delta \left(\phi \right)\left|\nabla \phi \right|dx+\alpha {\int }_{\Omega }\mathit{gH}\left(-\phi \right)dx$$(5) to generate a new external energy functional equation: (7) $${\epsilon }_{\mathit{unified}}\left(\phi \right)=\lambda {\int }_{\Omega }g\delta \left(\phi \right)\left|\nabla \phi \right|dx+\alpha {\int }_{\Omega }g\cdot \mathit{SPF}\left(I\left(x\right)\right)H\left(-\phi \right)dx.$$(7)

Consequently, the energy functional of the proposed modified level set can be represented as: (8) $$\epsilon \left(\phi \right)=\mu D_p\left(\phi \right)+{\epsilon }_{\mathit{unified}}\left(\phi \right).$$(8)

2.5. Evolution

Assume I₀ denotes the initial contour binary image of the live. Due to the existence of the dual-well function p(|▽φ|) in the DRLSE framework, by using binary step function, the level set the function G₀ can be initialized directly, and thus simplify the initialization of SPF. The initialization equation is as follows: (9) $$G_0=\mathit{\omega }\cdot \left(I_0-\frac{1}{2}\right)$$(9) where ω is a variable parameter for controlling the width of SDF. Since the image domain is discretely gridded, for the symbol distance function, there is at least one grid point on both sides of the zero level set contour. Therefore, the width of the symbol distance function is at least two [16].

Through the evolution of the initialized level set, the final refinement of liver segmentation is then achieved.

3. Experiment

In the training phase, we use 20 liver contrast-enhanced CT training datasets with corresponding segmentations, which are from MICCAI 2007 (www.sliver07.org). For the training dataset, the pixel size is: 0.55 × 0.55 − 0.80 × 0.80 mm, and the layer spacing is: 1.0-3.0 mm.

In the testing phase, 15 test datasets were used to validate the accuracy of the presented method. The test dataset and its corresponding “gold standard” segmentation results are from 3 D-IRCADb (https://www.ircad.fr/research/3dircadb/). For the test dataset, the pixel size is: 0.56 × 0.56 − 0.87 × 0.87 mm, with layer spacing: 1.0-4.0 mm. Each slice in all datasets contains flat pixels of 512 × 512.

In the experiments, two comparative studies are mainly implemented. One is the comparison on segmentation errors using the proposed patient-specific PA and standard PA, respectively, and another is a comparison on segmentation errors between our proposed method and two other methods.

The experiment is implemented with MATLAB R2010a, Intel i7-4770, 16GB RAM, and in order to evaluate the accuracy, 5 metrics were chosen in this research: VOE, RVD, ASD, RMSD, and MSD [1].

4. Result

According to the above five metrics, Table 1 gives the evaluation details by our proposed method on 15 3Dircadb datasets, in which the golden standard was provided by radiologist.

Table 1. Segmentation errors of our method by five metrics.

Subject	VOE (%)	RVD (%)	ASD (mm)	MSD (mm)	RMSD (mm)
1	6.47	3.51	1.23	18.14	2.15
2	9.32	−3.45	1.61	19.77	2.72
3	8.65	4.26	1.85	19.62	3.23
4	8.21	3.38	1.65	19.81	2.85
5	5.53	−0.85	1.14	14.13	1.96
6	7.76	−2.52	1.33	19.24	2.71
7	10.51	5.53	2.26	21.47	3.73
8	9.49	4.73	1.67	19.92	2.85
9	9.61	5.11	1.91	20.15	3.42
10	8.51	3.97	1.72	19.21	3.31
11	5.87	3.16	1.23	15.93	2.06
12	6.82	3.72	1.44	18.22	2.27
13	6.16	−1.79	1.08	11.39	1.71
14	7.35	4.35	1.67	18.83	2.72
15	8.97	4.51	1.46	18.53	2.61
Mean	7.95	2.51	1.54	18.29	2.69
Std. Dev.	1.48	2.91	0.31	2.52	0.56

To prove the superiority of patient-specific probabilistic atlas compared to standard atlas in our research, a comparative experiment is also implemented on the 15 3Dircadb datasets. The result is shown in Table 2. In addition, a comparison between the proposed method and two other approaches is shown in Table 3.

Table 2. Comparisons between standard atlas and our patient-specific PA on 15 3Dircadb datasets.

Method	VOE (%)	RVD (%)	ASD (mm)	MSD (mm)	RMSD (mm)
Standard Atlas	8.21  ±  1.72	2.75 ± 3.13	1.67 ± 0.54	21.11 ± 4.60	2.83 ± 0.47
Our PA	7.95  ±  1.48	2.51  ±  2.91	1.54  ±  0.31	18.29  ±  2.52	2.69  ±  0.56

Table 3. Comparisons with other methods with 3Dircadb datasets.

Method	VOE (%)	RVD (%)	ASD (mm)	MSD (mm)	RMSD (mm)
Li et al. [19]	9.15 ± 1.44^*	−0.07 ± 3.64	1.55 ± 0.39^*	28.22 ± 8.31	3.15 ± 0.98
Erdt et al. [20]	10.34 ± 3.11^*	1.55 ± 6.49	1.74 ± 0.59^*	26.83 ± 8.87	3.51 ± 1.16
Ours	7.95  ±  1.48	2.51  ±  2.91	1.54  ±  0.31	18.29  ±  2.52	2.69  ±  0.56

The statistically significant difference at a level of 0.05 is indicated by^*.

Figure 4 shows partial result of the proposed method (dark thick line) and manual segmentation (light thin line) on some representative slices. The black arrows indicated the liver tumor area. Comparing the results of our proposed method with ground truth, we can find that, even for slices containing tumors, our proposed method can still achieve accurate segmentation, with strong robustness.

Figure 4. Liver segmentation results of our proposed method.

For a same segmentation, the comparative 3D visualizations using our proposed method and ground truth were depicted in Figure 5(a,b), respectively. Based on the average ASD error, a 3 D visualizations of error between our result and ground truth result is also provided in Figure 5(c), in which, the ASD error was 1.41 mm.

Figure 5. 3D visualization of ground truth and our proposed method. (a) Ground truth. (b) Our segmentation. (c) Error image.

5. Discussion

It can be seen from Table 1 that, patient-specific PA obtained a higher accuracy on the five metrics. The main reason is that, hepatic diseases may affected some liver tissues and change its shape, meanwhile possible shape variations may also be lost while constructing such an atlas by standard atlas. To overcome these problems, standard atlas often requires a large-scale training datasets [21], while in this paper the weighted patient-specific probabilistic atlas is proposed based on the similarity between target testing image and multiple clusters, so that the obtained initial shape is much more closer to the shape and position of the testing image of the patient, by making full use of limited training samples.

From the Table 3, it can be seen that, our method outperforms the other two methods on VOE, ASD, RMSD and MSD, however, a worse result is obtained on RVD. The main reason for that is, although MLLR could retain possible liver contours to the greatest extent, it also causes over-segmentation, which inevitably results in a bigger RVD error. Moreover, an attempt was also made on liver image with tumor volume more than 50%, a poor MLLR was obtained by our proposed method, and resulting in a large segmentation error in the post fine-tune processing, which indicates the limitation of our proposed method when deployed on liver CT with larger tumors. In addition, the average running time for one scan was less than 7 minutes, which was mainly exhausted on solving the minimization problem of level set.

Furthermore, to evaluate whether the difference of segmentation accuracy between the proposed method and the other two methods was statistically significant, the paired t-test was conducted, which are based on two major metrics: VOE and ASD. As a result, our proposed method outperforms the other two technologies with statistically significant differences on VOE and ASD (p < 0.05).

6. Conclusion

In this paper, a Modified Distance Regularized Level Set-based liver segmentation method is proposed using patient-specific probabilistic atlas from abdominal CT. To some challenging clinical problems, the applicability of the proposed scheme is demonstrated, and results in good performance on both quantitative and qualitative analyses. Experiment shows that the proposed scheme can delineate liver boundaries that have a level similar to those obtained manually, and with good robustness.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported in part by the National Natural Science Foundation of China (No.61741106, 61702135, 61571158), the Natural Science Foundation of Heilongjiang Province (No. LH2019F203), the Fundamental Research Fundation for Universities of Heilongjiang Province (No. LGYC2018JQ004), the Innovative Personnel Project of Heilongjiang Province (No. UNPYSCT-2017089), the Higher Educational Science and Technology Project of Shandong Province (No. J18KB150), and the Innovative Personnel Project of Harbin (No. 2017RAQXJ205).